ベイズ統計モデリング入門

はじめに

統計モデリングとベイズ推論を行う上で、理解したい三つの概念

• モデル
• パラメータ
  
• データ

モデルを立て、与えられたデータを用いて、パラメータを推定し、未知のデータを予測したり現象の起こる原理を説明する。

データ

データはある種のランダムさを持って得られたものであると考える。

データの発生の仕方をうまく説明したい。 また推定されたパラメータに基づいて、将来的に得られるであろうデータを予測したい。

モデル

データの発生の仕方を説明するための大枠。 または分析する前提としての枠組み。

基本的にはモデルをどう設定するか、よいモデルを作れるかが分析者のやること。

パラメータ

モデルの細かな設定を決めるための数値などをいう。

確率

乱数やサンプリングという考え方について、Rやstanを用いながら親しみを持ってもらう。

例題

まずは上の統計モデリングの三つの概念について、特にパラメータの推定について理解するため、 次のような簡単な問題を考えよう。

手元にコインがある。これを投げて表が出るか、裏が出るかを予測したい。 すでに10回投げて7回表が出るということを観測している。

• モデルとして、コインはある一定の確率\(p\)で表が出て、\(1-p\)で裏が出るとする。
• この\(p\)がパラメータである。
• また10回投げて7回表が出た、というのが現在のデータである。

このデータを用いて\(p\)を推定する、というのがパラメータの推定。

推定されたパラメータを用いて、次に表が出るか裏が出るかを予測することができる。

どのようにパラメータ\(p\)を推定するか？

• 10回中7回なので\(p=0.7\)と推定
• コインだから表裏半々ぐらい出るはずなので、ちょっと表が出やすくて\(p=0.6\)と推定

など、色々考え方がある。

一つ目の推測は、最尤法という方法に基づいた結論と同じものになる。 最尤法の考え方を理解するために、次のような実験を行うことにしよう。

x <- rbinom(n = 1, size = 10, prob = 0.7)
x

## [1] 6

y <- rbinom(n = 1, size = 10, prob = 0.6)
y

## [1] 7

z <- rbinom(n = 1, size = 10, prob = 0.5)
z

## [1] 8

x1 <- rbinom(n = 100, size = 10, prob = 0.8)
x1

##   [1]  9  6  8 10  7  7  7  8  8 10  5  6  8  9  5  8  9  8  8 10  9  8  8
##  [24]  9  7  8  9 10  8  8  7  7  7  6  8  9  8  8  8  8  9  7  7  7  7  7
##  [47]  7  9  9  8 10  8  8  8  9  9  8  8  9  7  8  3  8  7  9  8  8  8  7
##  [70]  7  8  7  8  8  9 10  8  9  9  6  8  6  8 10  8  9 10  8  6  8  9  9
##  [93]  9  7  8  5  8  9  7  6

x2 <- rbinom(n = 100, size = 10, prob = 0.7)
x2

##   [1]  4  8  9  6  6  8  7  6  7  6  9  8  7  6  9  8  7  7  7  7  7  7  8
##  [24]  8  7  6  9  7 10  4  6  8  5  7  4  5  6  9  7  7  8  6  8  8  5  5
##  [47]  7  7  8  6  7  7  7  8  6  7  7  7  7  8  7  6  5  8  7  6  5  7  6
##  [70]  6  6  6  7  7  7  7  9  5  8  8  8  7  9  6  9  7  4  9  7  6  7  4
##  [93]  5  8  8  7  7 10  3  7

x3 <- rbinom(n = 100, size = 10, prob = 0.6)
x3

##   [1]  5  6  6  6  3  6  8  7  2  7  4  6  7  7  7  8  8  5  7  6  6  6  3
##  [24]  6  7  4  9  5  5  7  7  5  5  4  5  7  5  4  8  6  7  9  8  7  7  6
##  [47]  8  7  6  5  7  7  7  6  8  4  7  7  7  3  6  4  4  5  5  7 10  3  4
##  [70]  7  6  5  7  2  9  6  6  6  7  6  5  6  6  4  5  5  8  4  6  2  8  6
##  [93]  8  6  3  6  6  9  8  6

x4 <- rbinom(n = 100, size = 10, prob = 0.5)
x4

##   [1] 4 4 3 4 4 5 8 5 4 4 8 5 7 6 3 4 8 4 4 4 5 3 6 5 6 3 4 4 2 6 5 0 4 7 5
##  [36] 5 3 3 5 4 5 7 4 7 4 5 4 7 6 5 4 5 3 5 6 3 5 5 5 5 8 7 2 3 4 5 7 4 4 4
##  [71] 6 6 3 6 4 8 6 6 4 6 4 3 7 5 6 7 3 7 3 5 9 5 7 5 5 5 5 4 5 7

par(mfrow=c(2,2))
hist(x1, breaks = 0:10)
hist(x2, breaks = 0:10)
hist(x3, breaks = 0:10)
hist(x4, breaks = 0:10)

このようにパラメータ\(p\)の値によって、表の出る回数にばらつきがある。 この表のでる回数の割合をパラメータ\(p\)に対する尤度という。

今回のデータでは、表の出た回数が\(7\)回なので、 パラメータ\(p\)ごとに表が\(7\)回出る確率を計算し、それを\(L(p)\)と表す。 これを尤度関数という。

今回のデータに基づいて尤度関数を計算すると \[ L(p)=120~p^7(1-p)^3 \] となる。

最尤法では尤度関数が最大となる\(p\)をパラメータの推定値とする。 これのグラフをかくと、

f <- function(p){120*p^7*(1-p)^3}
plot(f, 0, 1)

であり、\(p=0.7\)の時に最大となることがわかる。

今回はパラメータ\(p\)を持つ二項分布モデルに基づいて、現象を説明している。

ではベイズ推論の枠組みでは、上の問題に対してどのような回答をするか？

パラメータ\(p\)をある決まった値として推定するのではなく、 ある程度幅を持ったなぜ分布として推定する。

分布にすると何が違うかを説明する。 例えば与えられたデータが

• 2回投げて1回表
• 100回投げて50回表

というふた通りの可能性を考えよう。 この時、最尤法でパラメータを推定するとどちらも\(p=0.5\)となる。 しかし、後者の方が信頼できそうだろう。

このような信頼度の違いを分布の形状として表現することができる。

例えば\(p\)の分布が以下の二つでどのように違うか検討せよ。

f1 <- function(p){choose(2,1) * p^2 * (1-p)^2}
f2 <- function(p){choose(100,50) * p^50 * (1-p)^50}

par(mfrow=c(1,2))
plot(f1, 0, 1)
plot(f2, 0, 1)

これがパラメータに分布を持たせることの一つの利点である。

ベイズ推論の考え方

パラメータはある決まった真の値があるのではなく、 何らかの確率分布にしたがっているとして推測を行う。

• 尤度関数
• 事前分布
• 事後分布

の三つの概念について理解する。

尤度関数

上ですでに出てきたように、パラメータを入力すると現在観測しているデータが発生する確率を計算する関数。

モデルとデータが与えられると関数が定まる。

事前分布と事後分布

いずれもパラメータの確率分布。 データを観測する前の分布のことを事前分布、データを観測した後の分布を事後分布という。

観測を通して事前分布を事後分布に取り替えることを、パラメータ更新、ベイズ更新などという。

パラメータの分布

ベイズ推論ではデータの発生を確率現象と捉えるのみでなく、パラメータも確率的な現象と考える。

確率分布とベイズの定理

ここからは、ベイズ推論の枠組みを理解するために以下の設定を通して、 確率分布、特にパラメータ分布の更新について理解しよう。

設定

外からは見分けのつかない二つの箱A, Bの中に赤いボールと白いボールがいくつか入っている。

以下の手順に従ってボールを取り出し赤いボールなら賞金1万円、白いボールなら罰金1万円もらえるというゲームを考えよう。

• まずはじめにどちらかの箱を選ぶ。

• 次に観察期間として選んだ箱から規定の回数、ボールを取り出し色を確認して戻すという操作を行う。

• そののち、実際にゲームに挑戦するか決める。

• 挑戦する場合は先ほど選んだ箱の中からボールを一つとり、色に応じて賞金を得る。

さて、どのような戦略を立てるべきか。

例題1

• 箱Aには赤いボールが5個、白いボールが0個
• 箱Bには赤いボールが0個、白いボールが5個

入っているということがあらかじめわかっているとしよう。

この場合には、観察期間として1回ボールを確認すれば、選んだ箱がAなのかBなのか完全に決定することができる。

従ってそれに応じて最終的に挑戦するかどうか決定すればよい。

例題2

• 箱Aには赤いボールが3個、白いボールが1個
• 箱Bには赤いボールが1個、白いボールが3個

入っていることがわかっているとしよう。

また観察期間として、一度だけボールの色を確認できるとする。

Rの乱数を利用して実験してみよう。

rbinom(n = 10, size = 1, prob = 0.75)

##  [1] 0 0 1 1 0 0 1 1 1 1

rbinom(n = 10, size = 1, prob = 0.25)

##  [1] 1 1 1 0 0 1 1 1 0 0

赤いボールが出た場合には箱Aと推定し、白いボールが出た場合には箱Bと推定するのがもっともらしいであろう。

さらにそれを踏まえた上で、最終的に挑戦するかを決定できる。

このような判断を数学的に考えていこう。 まずは

• 確率分布
• 条件付き確率
• 同時確率
• 周辺確率
• 尤度関数
• ベイズの定理

という概念を順番に理解していこう。

確率分布

\(X\)という出来事が起こる確率を\(P(X)\)とかく。 例えば \[ P(\mbox{コインを投げて表が出る})=\dfrac{1}{2} \] のように書く。 これを省略して \[ P(\mbox{表})=\dfrac{1}{2} \] などとも書く。

今回の問題では箱がAなのかBなのか、ボールが赤なのか白なのか、といった出来事について確率を考える。 \[ P(A)=\dfrac{1}{2} \] \[ P(B)=\dfrac{1}{2} \] と設定しよう。 二つの箱からランダムに選ぶのでこのように設定するのが妥当だろう。

確率変数と確率分布

上のようにある状況で起こりうる出来事全体が\(X_0, X_1, \ldots, X_n\)である時、 \(P(X_1)\)から\(P(X_n)\)にたいして0以上1以下の数を与える。

このようなものを（離散的な）確率分布という。 この時、合計が1になるように決めなければならない。

例1

「コインを1回投げる」という操作に注目すると、起こりうる結果は表が出るか裏が出るかどちらか。 この結果を表す確率変数\(X\)を導入する。

\(P(\mbox{表})\)と\(P(\mbox{裏})\)を決めればこの現象を記述できる、と考える。 この二つの値が0以上1以下で、合計が1になるように配分してやる。 例えば \[ P(X=\mbox{表})=\dfrac{1}{2},~P(X=\mbox{裏})=\dfrac{1}{2} \] としてもよいし、 \[ P(X=\mbox{表})=\dfrac{1}{3},~P(X=\mbox{裏})=\dfrac{2}{3} \] としてもよい。 極端な場合には \[ P(X=\mbox{表})=1,~P(X=\mbox{裏})=0 \] としてもよい。

例2

「サイコロを1回投げる」という操作に注目すると、起こりうる結果は1,2,3,4,5,6のどれかが出るということ。 サイコロの目を表す確率変数を\(X\)とする。

この場合は \[ P(X=1),~P(X=2),~P(X=3),~P(X=4),~P(X=5),~P(X=6) \] を決めてやる。 例えば \[ P(X=1)=\dfrac{1}{6},~P(X=2)=\dfrac{1}{6},~P(X=3)=\dfrac{1}{6},~P(X=4)=\dfrac{1}{6},~P(X=5)=\dfrac{1}{6},~P(X=6)=\dfrac{1}{6} \] としてもよいし、 \[ P(X=1)=\dfrac{1}{10},~P(X=2)=\dfrac{2}{10},~P(X=3)=\dfrac{3}{10},~P(X=4)=\dfrac{1}{10},~P(X=5)=\dfrac{2}{10},~P(X=6)=\dfrac{1}{10} \] としてもよい。

条件付き確率と同時分布

次にボールを一つ取り出した時、赤であるか白であるかの確率を計算したいが、このままでは難しい。

なぜならどちらの箱を選んだかによって赤いボールの出やすさが異なるからである。

このように箱の種類とボールの色といういくつかの条件を組み合わせる場合には、 条件付き確率や同時分布を用いて記述する。

例えば選んだ箱がAだという条件のもと、赤いボールを取り出す確率を表現したい。

\(Y\)という条件のもとでの\(X\)の確率を \[ P(X\mid Y) \] と表し、これを条件付き確率という。

箱の種類を表す確率変数を\(\theta\)とし、取り出したボールの色を表す変数を\(X\)とする。

Aの箱という条件のもとで赤いボールを取り出す確率は \[ P(X=\mbox{赤}\mid \theta=A)=\frac{3}{4} \] であり、Aの箱という条件のもとで白いボールを取り出す確率は \[ P(X=\mbox{白}\mid \theta=A)=\frac{1}{4} \] となる。 この二つの合計は1になっている。 よってこれは確率分布を与えていることに注意しよう。

同様に箱Bを選んだという条件のもとでのボールの色に関する条件付き確率を表現すると、 \[ P(X=\mbox{赤}\mid \theta=B)=\frac{1}{4},~ P(X=\mbox{白}\mid \theta=B)=\frac{3}{4} \] と書くことができる。 これも確率分布を与えている。

同時分布と周辺分布

箱を選びボールを取り出す、という二つの条件をまとめて考えると、起こりうる出来事としては

• 箱Aから赤いボール
• 箱Aから白いボール
• 箱Bから赤いボール
• 箱Bから白いボール

の4パターンがある。 このように二つ以上の条件をまとめて考えた確率を同時確率という。

これらの確率を計算すると、 \[ P(X=\mbox{赤}\cap \theta=A)=\frac{3}{4}\times\frac{1}{2}=\frac{3}{8} \] \[ P(X=\mbox{白}\cap \theta=A)=\frac{1}{4}\times\frac{1}{2}=\frac{1}{8} \] \[ P(X=\mbox{赤}\cap \theta=B)=\frac{1}{4}\times\frac{1}{2}=\frac{1}{8} \] \[ P(X=\mbox{白}\cap \theta=B)=\frac{3}{4}\times\frac{1}{2}=\frac{3}{8} \] と計算できる。 4つを合計すると1になるので、これも確率分布である。

また\(P(\theta = A), P(\theta = B)\)や\(P(X= \mbox{赤}), P(X=\mbox{白})\)と行ったように、どちらかの条件のみを考えた確率分布を周辺分布という。 この場合には \[ P(X=\mbox{赤})=P(X=\mbox{赤}\cap \theta=A)+P(X=\mbox{赤}\cap \theta=B)=\frac{1}{2} \] \[ P(X=\mbox{白})=P(X=\mbox{白}\cap \theta=A)+P(X=\mbox{白}\cap \theta=B)=\frac{1}{2} \] として計算できる。 クロス集計表を思い浮かべるとよい。

ベイズの定理

クロス集計表をじっと見ていると、次のような関係が成り立つことがわかる。

\[ P(X\mid Y)P(Y)=P(X\cap Y)=P(Y\mid X)P(X) \]

これを式変形して得られる次の公式をベイズの定理という。 \[ P(X\mid Y)=\frac{P(Y\mid X)P(X)}{P(Y)} \] これを用いることで、\(P(A\mid \mbox{赤})\)と\(P(B\mid \mbox{赤})\)のいずれが大きいかを計算することができる。

言い換えると、赤いボールを一個とったあとで、目の前の箱がAであるのかBであるのか、どちらの確率が高いかを計算できる。

例題1の場合で言えば、 Aの箱を選んだという条件のもとで赤いボールが出る確率は \[ P(X=\mbox{赤}\mid \theta=A)=1 \] であり、Aの箱を選んだという条件のもとで白いボールが出る確率は \[ P(X=\mbox{白}\mid \theta=A)=0 \] となる。 またBの箱から取ったという条件のもとでの赤白のそれぞれの確率は \[ P(X=\mbox{赤}\mid \theta=B)=0,~P(X=\mbox{白}\mid \theta=B)=1 \] となる。

同時分布は \[ P(\theta=A\cap X=\mbox{白})=0,~P(\theta=A\cap X=\mbox{赤})=\frac{1}{2},~P(\theta=B\cap X=\mbox{白})=\frac{1}{2},~P(\theta=B\cap X=\mbox{赤})=0 \] となり、さらに \[ P(X=\mbox{赤})=P(\theta=A\cap X=\mbox{赤})+P(\theta=B\cap X=\mbox{赤})=\frac{1}{2} \] \[ P(X=\mbox{白})=P(\theta=A\cap X=\mbox{白})+P(\theta=B\cap X=\mbox{白})=\frac{1}{2} \] であることがわかる。

さて、赤いボールを取り出したとき、選んだ箱がどちらの箱の確率が高いだろうか？

\(P(\theta=A\mid X=\mbox{赤})\)と\(P(\theta=B\mid X=\mbox{赤})\)の大小を比較しよう。

ベイズの定理を使えば、 \[ P(\theta=A\mid\mbox{赤})=\frac{P(\mbox{赤}\mid \theta=A)P(\theta=A)}{P(\mbox{赤})}=1\times\frac{1}{2}\div\frac{1}{2}=1 \] \[ P(\theta=B\mid\mbox{赤})=\frac{P(\mbox{赤}\mid \theta=B)P(\theta=B)}{P(\mbox{赤})}=0\times\frac{1}{2}\div\frac{1}{2}=0 \]

と計算できる。

例題2の場合を考えよう。

1回ボールを取り出すと赤いボールがでた、という事実を観測データ\(X\)としよう。

すると、赤いボールが出る確率をパラメータの関数として表現すると以下のようになる。

\[ L(A)=P(X=\mbox{赤}\mid \theta=A)=\frac{3}{4} \] \[ L(B)=P(X=\mbox{赤}\mid \theta=B)=\frac{1}{4} \] この\(L\)が尤度関数である。

ベイズの定理を使えば、 \[ P(\theta=A\mid X=\mbox{赤})=\frac{P(X=\mbox{赤}\mid \theta=A)P(\theta=A)}{P(X=\mbox{赤})}=\frac{L(A)P(\theta=A)}{P(X=\mbox{赤})}=\frac{3}{4}\times\frac{1}{2}\div\frac{1}{2}=\frac{3}{4} \] \[ P(\theta=B\mid X=\mbox{赤})=\frac{P(X=\mbox{赤}\mid \theta=B)P(\theta=B)}{P(X=\mbox{赤})}=\frac{L(B)P(\theta=B)}{P(X=\mbox{赤})}=\frac{1}{4}\times\frac{1}{2}\div\frac{1}{2}=\frac{1}{4} \]

と計算できる。

これが事後分布。

この事後分布を用いて、次のボール\(Y\)が赤である確率を計算すると、 \[ P(Y=\mbox{赤})=P(Y=\mbox{赤}\mid \theta=A)P(\theta=A\mid X=\mbox{赤})+P(Y=\mbox{赤}\mid \theta=B)P(\theta=B\mid X=\mbox{赤})=\frac{3}{4}\times\frac{3}{4}+\frac{1}{4}\times\frac{1}{4}=\frac{10}{16} \] となる。

さてここまでの言葉遣いを用いて、ベイズ推測の枠組みを改めて説明すると、以下のようになる。

\(P(\theta=A), P(\theta=B)\)という確率分布を事前分布、\(P(\theta=A\mid\mbox{白}), P(\theta=B\mid\mbox{白})\)という確率分布を事後分布という。

• モデルは確率分布\(P(\mbox{赤}\mid \theta=A), P(\mbox{白}\mid \theta=A)\)という確率分布から決まるもの
• パラメータは箱の種類AまたはB
• データは観察期間にでたボールの記録

ということになる。

例題3

次に箱の設定は前と同じで観察できる回数が5回になったとしよう。

前回とはデータとして想定する出来事の起こり方が変わってくる。

• 例題2ではボールを取り出して赤が出るか白が出るか
• 例題3ではボールを5回取り出して何回赤が出るか

という点が異なる。

このような出来事の確率分布を記述するため

• 二項分布

について理解する。

rbinom(n = 1, size = 5, prob = 0.75)

## [1] 4

rbinom(n = 1, size = 5, prob = 0.25)

## [1] 2

今回の問題設定は、このような数を見たときにprobの値がどのようになっているか推測する、ということができる。

二項分布

表と裏がそれぞれ\(\dfrac{1}{2}\)の確率で出るコイン、つまり\(P(\mbox{表})=P(\mbox{裏})=\dfrac{1}{2}\)であるようなコインを\(N\)回投げる、 という操作を考えよう。 \(N=1\)のとき、表が出るか裏が出るかどちらかで、表が\(1\)回出る確率、\(0\)回出る確率はどちらも\(\dfrac{1}{2}\)である。

\(N=2\)のとき、表表、表裏、裏表、裏裏の4パターンあり、 \[ P(\mbox{表表})=P(\mbox{表裏})=P(\mbox{裏表})=P(\mbox{裏裏})=\frac{1}{2}\times\frac{1}{2}=\frac{1}{4} \] と全て等しい。これは1回目の確率と2回目の確率をかけて計算できる。

出方のパターンではなく、表の出る回数\(X\)の分布を考えよう。 \(N=2\)の時は表の出る回数は0,1,2のいずれかなので\(P(X=0),~P(X=1),~P(X=2)\)という確率分布を与える。 上で計算したことから \[ P(X=0)=\frac{1}{4},~P(X=1)=\frac{1}{4}+\frac{1}{4}=\frac{1}{2},~P(X=2)=\frac{1}{4} \] となることがわかる。 \(X=1\)となるのは2通りあって、それらの確率の和を考えればよい。

\(N=3\)ならば、表裏のでる出方は表表表、表表裏、表裏表、表裏裏、裏表表、裏表裏、裏裏表、裏裏裏、と八通りあり、いずれも確率\(\dfrac{1}{2}\)である。 従って \[ P(X=0)=\frac{1}{8},~P(X=1)=\frac{1}{8}+\frac{1}{8}+\frac{1}{8}=\frac{3}{8},~P(X=2)=\frac{1}{8}+\frac{1}{8}+\frac{1}{8}=\frac{3}{8},~P(X=3)=\frac{1}{8} \] と計算できる。

\(N=4\)の場合、表裏の出るパターンを全て列挙し、表の出る回数の確率分布を計算しよう。

一般に\(N\)回コインを投げた時、表が\(k\)回出るパターンは \[ {}_NC_k=\dfrac{N!}{k!(N-k)!} \] 通りある。

例えば \[ {}_2C_1=\frac{2}{1}=2,~_4C_2=\frac{4!}{2!\times2!}=6 \] などとなる。

したがって\(N\)回コインを投げた時の、表が出る回数の確率分布は \[ P(X=k)=~_NC_k\left(\frac{1}{2}\right)^k\left(\frac{1}{2}\right)^{N-k} \] と計算できる。

より一般に1回コインを投げた時に表の出る確率\(p\)とする。 つまり \[ P(\mbox{表})=p,~P(\mbox{裏})=1-p \] とする。 この時\(N\)回コインを投げた時の、表が出る回数の確率分布は \[ P(X=k)=~_NC_k~p^k(1-p)^{N-k} \] となる。 これが二項分布である。

例題4

箱の個数がA, Bで異なるとしよう。 例えばAの箱が3個、Bの箱が2個あるとする。

この時は、パラメータの事前分布を \[ P(\theta=A)=\dfrac{3}{5},~P(\theta=B)=\dfrac{2}{5} \] と設定する。

5回ボールを取り出した結果、赤いボールが4個、白いボールが1個、というデータが与えられたとしよう。

確率変数\(X\)を赤いボールの個数とする。

このとき尤度関数は、 \[ L(A)=P(X=4\mid \theta=A)={}_5C_4~\left(\frac{3}{4}\right)^4\left(\frac{1}{4}\right)^1=\frac{5\times3^4}{4^5} \] \[ L(B)=P(X=4\mid \theta=B)={}_5C_4~\left(\frac{3}{4}\right)^1\left(\frac{1}{4}\right)^4=\frac{5\times3}{4^5} \] であり、 \[ P(X=4)=P(X=4\mid \theta=A)P(\theta=A)+P(X=4\mid \theta=B)P(\theta=B)=\frac{5\times3^4}{4^5}\times\frac{3}{5}+\frac{5\times3}{4^5}\times\frac{2}{5}=\frac{3^5+6}{4^5} \] となる。

したがってベイズの定理を使えば、事後分布は \[ P(\theta=A\mid X=4)=\frac{P(X=4\mid \theta=A)P(\theta=A)}{P(X=4)}=\frac{L(A)P(\theta=A)}{P(X=4)}=\frac{5\times3^4}{4^5}\times\frac{3}{5}\div\frac{3^5+6}{4^5}=\frac{3^5}{3^5+6} \] \[ P(\theta=B\mid X=4)=\frac{P(X=4\mid \theta=B)P(\theta=B)}{P(X=4)}=\frac{L(B)P(\theta=B)}{P(X=4)}=\frac{5\times3}{4^5}\times\frac{2}{5}\div\frac{3^5+6}{4^5}=\frac{6}{3^5+6} \]

と計算できる。

この事後分布を用いて、次のボール\(Y\)が赤である確率を計算すると、 \[ P(Y=\mbox{赤})=P(Y=\mbox{赤}\mid \theta=A)P(\theta=A\mid X=4)+P(Y=\mbox{赤}\mid \theta=B)P(\theta=B\mid X=4)=\frac{3}{4}\times\frac{3^5}{3^5+6}+\frac{1}{4}\times\frac{6}{3^5+6} \] となる。

例題5

次に箱の種類がA, B, C, Dと4種類あり、それぞれ赤白のボールの内訳が、

• 箱Aは赤いボールが1個、白いボールが4個
• 箱Bは赤いボールが2個、白いボールが3個
• 箱Cは赤いボールが3個、白いボールが2個
• 箱Dは赤いボールが4個、白いボールが1個

となっているとしよう。 またこの箱が全て一つずつあるとする。

パラメータの事前分布としては \[ P(\theta= A)=P(\theta=B)=P(\theta= C)=P(\theta= D)=\frac{1}{4} \] とすればよい。

5個ボールを取り出して赤いボールを3個白いボールを2個取り出した、というデータが与えられたとしよう。

確率変数\(X\)を赤いボールの個数とする。

このとき尤度関数は、 \[ L(A)=P(X=3\mid \theta= A)={}_5C_3~\left(\frac{1}{5}\right)^3\left(\frac{4}{5}\right)^2 \] \[ L(B)=P(X=3\mid \theta=B)={}_5C_3~\left(\frac{2}{5}\right)^3\left(\frac{3}{5}\right)^2 \] \[ L(C)=P(X=3\mid \theta= C)={}_5C_3~\left(\frac{3}{5}\right)^3\left(\frac{2}{5}\right)^2 \] \[ L(D)=P(X=3\mid \theta= D)={}_5C_3~\left(\frac{4}{5}\right)^3\left(\frac{1}{5}\right)^2 \] であり、 \[ P(X=3)=P(X=3\mid \theta = A)P(\theta= A)+P(X=3\mid \theta=B)P(\theta=B)+P(X=3\mid \theta= C)P(\theta= C)+P(X=3\mid \theta = D)P(\theta = D) \]

\[ =\frac{1}{4}\times10\times\frac{1}{5^5}\times(4^2+2^3\times3^2+3^3\times2^2+4^3)=\frac{26}{5^3} \] となる。

したがってベイズの定理を使えば事後分布を、 \[ P(\theta= A\mid X=3)=\frac{P(X=3\mid \theta = A)P(\theta= A)}{P(X=3)}=\frac{L(A)P(\theta= A)}{P(X=3)}=10\times\frac{4^2}{5^5}\times\frac{1}{4}\div\frac{26}{5^3}=\frac{4}{65} \] \[ P(\theta=B\mid X=3)=\frac{P(X=3\mid \theta=B)P(\theta=B)}{P(X=3)}=\frac{L(B)P(\theta=B)}{P(X=3)}=10\times\frac{2^3\times3^2}{5^5}\times\frac{1}{4}\div\frac{26}{5^3}=\frac{18}{65} \] \[ P(\theta= \theta= C\mid X=3)=\frac{P(X=3\mid \theta= C)P(\theta= C)}{P(X=3)}=\frac{L(C)P(\theta= C)}{P(X=3)}=10\times\frac{3^3\times2^2}{5^5}\times\frac{1}{4}\div\frac{26}{5^3}=\frac{27}{65} \] \[ P(\theta= D\mid X=3)=\frac{P(X=3\mid \theta= D)P(\theta= D)}{P(X=3)}=\frac{L(D)P(\theta= D)}{P(X=3)}=10\times\frac{4^3}{5^5}\times\frac{1}{4}\div\frac{26}{5^3}=\frac{16}{65} \]

と計算できる。

この事後分布を用いて、次のボール\(Y\)が赤である確率は \[ P(Y = \mbox{赤})=P(Y= \mbox{赤}\mid \theta= A)P(\theta= A\mid X=3)+P(Y= \mbox{赤}\mid \theta=B)P(\theta= B\mid X=3)+P(Y= \mbox{赤}\mid \theta=C)P(\theta= C\mid X=3)+P(Y= \mbox{赤}\mid \theta=D)P(\theta= D\mid X=3)= \] により計算できる。

例題6

箱AからDの中身は前と同じだが、 今回は箱の個数がどのようになっているかわからないとする。

特に何も情報がない場合、全ての箱の割合が均等である一様分布を設定するべきであろう。

また、事前に情報がある場合や、例えば賞金や罰金の額が均等でないといった場合、何らかの情報を持った事前分布を設定することもありうる。

このようにパラメータの事前分布と呼ばれる確率分布を自分で設定する必要がある。 多くの場合、事前分布は一様分布を設定する。 これを無情報事前分布といったりする。

例題7

さらに箱の種類が101種類、型番0番から100番まであり、赤いボールの割合\(p\)が型番と一致するとする。

さらに、各番号の箱の個数が一定であるかはわからない状況であるとする。

この場合、モデルが二項分布であることは同じであるが、事前分布の設定を変えることになる。

まず尤度関数を計算しよう。

型番\(p\)の箱を選んでいるとしたとき、\(n\)回ボールを取り出したうち\(k\)回赤いボールになる確率は \[ L(p)=P(X=k\mid 箱~p)={}_nC_k~p^k(1-p)^{n-k} \] と計算できる。

例えば、5回中4回赤いボールというデータが与えられたとすると、 \[ L(p)=P(X=4\mid 箱~p)={}_5C_4~p^4(1-p)^1 \]

さて、事前分布は各\(p\)ごとに\(p\)番の箱がいくつあるか？を事前に設定しておき、 各\(p\)ごとに値を定めておくものだった。

今回の設定では箱の種類が多いので、連続的な分布を用いて近似的に考えることにしよう。

連続的な分布の場合は、\(p\)の取りうる値が無数にあるので、 各\(p\)の箱の割合を与えるのではなく、\(p\)付近の密度を与える関数を考える。 それを確率密度関数という。

例えば事前になんの情報もない場合には\(p=0\)から\(p=1\)までが均等にあることを表す一様分布を事前分布として用いる。 一様分布の確率密度関数は \[ f(p)=1~(0\leq p \leq1) \] である。

上と同じようにベイズの定理を用いて事後分布を計算すると、 \[ f(p\mid X=4)=\frac{L(p)f(p)}{P(X=4)}=\frac{5}{P(X=4)}p^4(1-p)^1 \] となる。

ここで\(P(X=4)\)を計算するためには積分の計算が必要になる。

例題8

箱に関する設定は上の例題と同じとしよう。

なんらかの情報がある場合には事前分布を自由に設定できる。

ここではベータ分布と呼ばれる確率分布を使うことにしよう。 ベータ分布の確率密度関数は \[ f(p)=\frac{1}{C}~p^{\alpha-1}(1-p)^{\beta-1} \] である。 ここで\(\alpha, \beta\)はベータ分布のパラメータである。

また分母の\(C\)は面積を1にするため（確率分布の条件）の定数なので、今回はあまり気にしないでよい。

この\(x\)の関数を密度関数に持つ確率分布をベータ分布という。 例えば\(Be(1,1)\)は一様分布である。

ベータ分布からサンプリングしてみよう。 またalpha, betaを変えると分布の様子がどのように変わるか観察してみよう。

n <- 10000
alpha <- 3
beta <- 4
x<-rbeta(n = n, shape1 = alpha, shape2 = beta)
f <- function(x){x^(alpha-1)*(1-x)^(beta-1)*n*0.05/beta(alpha, beta)}
hist(x)
plot(f,0,1, add=TRUE)

事前分布をベータ分布とした時、 事後分布は \[ f(p\mid X=k)=L(p)\times f(p)=p^{\alpha-1}(1-p)^{\beta-1}\times p^k(1-p)^{N-k}=p^{\alpha+k-1}(1-p)^{\beta+N-k-1} \] となるので、またベータ分布になる。

連続型確率分布

改めて連続型確率分布の考え方について整理しよう。

データやパラメータなどの値が有限でない、あるいはそのように近似できる場合を考える。 例えば、色々なものの大きさや量を測る場合を考えよう。 このようなものを連続型の確率変数という。

それに対して数え上げられるだけの可能性しかとらないものを離散型の確率変数という。

上の例題で言えば箱の種類であったり、ある出来事が何回起こるかといった整数の値で表されるものなどである。

離散型の確率変数の場合、確率分布とは起こりうる結果が\(X=x_1,x_2,\ldots,x_n\)に対して、 \[ P(X=x_1),~P(X=x_2),\ldots,P(X=x_n) \] に\(0\)以上の数を割り当て、それらの和が\(1\)となるように分布させたものだった。

一方で連続型の確率変数の場合、飛び飛びの値ではなくぎっしり詰まった（ある範囲の）数すべてを取ることになる。 したがって、ある決まった値をとる確率\(P(X=k)\)は\(0\)でないとおかしなことになってしまう。

連続分布の場合、\(P(a\leq X\leq b)\)のように範囲を指定して確率を計算する。

一般に、確率変数\(X\)が連続分布に従うとき、密度関数と呼ばれる関数\(f(x)\)があって\(P(a\leq X\leq b)\)は\(f(x)\)と\(x\)軸で囲まれた部分の面積として与えられる。 この面積を表す記号として

\[ \int^b_af(x)dx \] を定める。

例えばあとで見る標準正規分布の場合、 その密度関数は \[ f(x)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2}) \] で与えられ、そのグラフは

mu <- 0
sigma <- 1
f <- function(x){(1/sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2))}
plot(f,-4,4)

ここから乱数をサンプリングし、ヒストグラムをかくと、

N <- 100000
x <- rnorm(N)
par(mfrow=c(2,2))
hist(x,breaks = (0:20)/2-5)
hist(x,breaks = (0:40)/4-5)
hist(x,breaks = (0:80)/8-5)
hist(x,breaks = (0:160)/16-5)

となる。 幅を細くしていけば、密度関数の形に近づいていくことが見える。

ヒストグラムでは高さが\(x\)がある範囲にあるサンプルの個数で、全体の個数に対する割合が確率であると考えられる。

このように見れば、確率密度関数とその面積を用いて確率を表現できることが納得できる。

連続型確率分布に対するベイズの定理

ベイズ推論に置いて扱う事象は、

• パラメータがどのような条件をみたすか
• 得られたデータがどのようなものであるか

の大きく二つがある。

例えば先ほどまでの例題で言えば

• 箱の名前、種類を表す変数\(X\)
• 赤いボールを何個取り出したかを表す変数\(\theta\)

の二つがあった。

このように二つの連続型確率変数\(X, Y\)を同時に考えたい場合、 二変数の密度関数\(f(x,y)\)を与えることで確率分布を記述する。このとき \[ P(a\leq X\leq b, c\leq Y\leq d)=\int^b_a\int^d_cf(x,y)dxdy \]

として確率分布が与えられる。

\[ f(x)=\int^\infty_{-\infty}f(x,y)dy \]

とすれば\(X\)の密度関数が得られる。 これは離散分布の場合に

\[ P(X=x)=\sum_iP(X=x, Y=y_i) \]

としたのと同じこと。

またこの場合に条件付き分布は密度関数

\[ f(y\mid x)=\frac{f(x,y)}{f(x)} \]

により与えられる。

連続分布に対してもベイズの定理が成り立つ。 \[ f(y\mid x)=\frac{f(x\mid y)f(y)}{f(x)} \]

連続型確率分布の例

以下ではいくつかの連続型確率分布の例とそれを用いたベイズ推定の例題を紹介する。

• 一様分布
• 正規分布
• ガンマ分布

一様分布

ある範囲の実数が一様に現れるような確率分布を一様分布という。

例えば\(0\)以上\(1\)以下の実数が一様に現れるような確率分布であれば、その確率密度関数は \[ f(x)= \begin{cases} 1 & 0\leq x\leq 1\\ 0 & x<0, 1<x \end{cases} \] となる。 このとき、例えば\(P(0.5\leq x\leq 0.6)=0.1\)と計算できる。

正規分布

正規分布は次の密度関数\(f(x)\)で定まる確率分布。 \[ f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\mu)^2}{2\sigma^2}) \]

ここで\(\mu\)と\(\sigma\)は正規分布のパラメータ。

確率密度関数のグラフは\(x=\mu\)がピークで左右対称。 この幅の広がり具合を決めるのが\(\sigma\)であり、

パラメータ\(\mu, \sigma\)の正規分布に従う確率変数は

• 平均\(\mu\)
• 分散\(\sigma^2\)

である。

正規分布は誤差の分布として使われる。 色々と良い性質を持ち、もっとも基本的な連続分布。

mu <- 0
sigma <- 1
f <- function(x){(1/sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2))}
plot(f,-3,3)

確率はRで計算できる

pnorm(1)

## [1] 0.8413447

正規分布からのサンプリングは

mu <- 0
sigma <- 1
N <- 1000
x <- rnorm(n = N, mean = mu, sd = sigma)
hist(x)

ガンマ分布

ガンマ分布は次のような密度関数を持つ。 \[ f(x)=\frac{\theta^k}{\Gamma(k)}e^{-x\theta}x^{k-1} \]

theta <- 1
k <- 5
f <- function(x){exp(-x*theta)*x^{k-1}*theta^k/gamma(k)}
plot(f,0,10)

ガンマ分布のパラメータは\(k\)と\(\theta\)である。 ガンマ分布に従う確率変数は

• 平均は\(k/\theta\)
• 分散は\(k/\theta^2\)

となる。

連続型確率分布を用いたベイズ推定

それでは改めて連続型確率分布を用いたベイズ推定の問題を考えてみよう。

例題：ポアソン分布のパラメータ推定

ある地域での一日の事故件数\(X\)を予測したい。

これをポアソン分布モデルを用いて予測してみよう。 データとして、これまで10日間の事故件数データが、

x<-c(8,4,5,5,7,9,7,5,5,4)

と与えられているとしよう。

ポアソン分布はパラメータ\(\lambda\)をもつ離散型確率分布である。 パラメータ\(\lambda\)のポアソン分布に従う確率変数\(X\)は

\[ P(X=k)=e^{-\lambda}\frac{\lambda^k}{k!} \] と確率を計算でき、

• 平均\(\lambda\)
• 分散\(\lambda\)

である。

ポアソン分布からサンプリングした変数のヒストグラムを、\(\lambda\)の値を変えて観察してみよう。

N <- 10000
lambda <- 3
x <- rpois(n = N, lambda = lambda)
hist(x,breaks=0:max(x))

今回のデータから尤度関数を計算すると、 \[ L(\lambda) =e^{-\lambda}\frac{\lambda^{x_1}}{x_1!} \times\cdots\times e^{-\lambda}\frac{\lambda^{x_n}}{x_n!} =\prod_{i=1}^Ne^{-\lambda}\frac{\lambda^{x_i}}{x_i!} =e^{-10\lambda}\frac{\lambda^{59}}{x_1!\cdots x_n!} \]

Rで尤度関数を計算すると、

x<-c(8,4,5,5,7,9,7,5,5,4)
L <- function(lambda){
  lik <- 1
  for (i in 1:length(x)){
    lik <- lik * exp(-lambda)*lambda^x[i]/factorial(x[i])
  }
  return(lik)
}

パラメータ\(\lambda\)の事前分布として、ガンマ分布を用いることにする。

事前分布をパラメータ\(k=5, \theta=1\)のガンマ分布としよう。 \[ f(\lambda)=\frac{\theta^5}{\Gamma(5)}e^{-\lambda}\lambda^{4} \]

すると、事後分布は \[ e^{-11\lambda}\lambda^{63} \] により決まるガンマ分布になる。 つまりパラメータ\(k=64, \theta=11\)のガンマ分布である。

\[ p(\lambda\mid D)=\frac{9^{64}}{64!}e^{-9\lambda}\lambda^{63} \]

この事後分布を用いて、一日の事故件数を予測することができる。 \[ P(X=x)=\int_{-\infty}^{\infty} e^{-\lambda}\frac{\lambda^x}{x!}p(\lambda\mid D)d\lambda \]

例題

上の例題ではパラメータの分布に連続型確率分布を用いたが、データの分布にも連続型分布を用いることができる。

そのような例題として、正規分布に従うデータのベイズ推定の例題を用意した。

正規分布の例題

工業製品の規格を調査する。 ある工場で作られた製品の大きさが、以下のようであった。

x <- c(102.3,100.4,99.3,100.2,100.3,99.4,101.5,99.3,98.9,100.1)

製品が規格通り製造されているか調査したい。

正規分布モデルを用いる。 パラメータは平均\(\mu\)と標準偏差\(\sigma\)で、これに対して上のデータから尤度関数を計算すると \[ L(\mu,\sigma) =\frac{1}{\sqrt{2\pi}\sigma}\exp(\frac{-(x_1-\mu)^2}{2\sigma^2}) \times\cdots\times \frac{1}{\sqrt{2\pi}\sigma}\exp(\frac{-(x_N-\mu)^2}{2\sigma^2}) =\prod_{i=1}^N\frac{1}{\sqrt{2\pi}\sigma}\exp(\frac{-(x_i-\mu)^2}{2\sigma^2}) \] である。

x <- c(102.3,100.4,99.3,100.2,100.3,99.4,101.5,99.3,98.9,100.1)
L <- function(mu, sigma){
  L <- 1
  for(i in 1:10){
    L <- L*dnorm(x[i],mu,sigma)
  }
}

さてここでは

• 平均の事前分布として、平均100, 分散1の正規分布
• 分散の事前分布として、平均1, 分散1の逆ガンマ分布

を用いることにしよう。

\[ p_1(\mu)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(\mu-100)^2}{2}) \]

\[ p_2(\sigma^2)=4(\sigma^2)^{-4}\exp(-\frac{2}{\sigma^2}) \]

すると事後分布は \[ p(\mu,\sigma^2\mid D)=L(\mu,\sigma^2)p_1(\mu)p_2(\sigma^2) \] \[ =\prod_{i=1}^N\frac{1}{\sqrt{2\pi}\sigma}\exp(\frac{-(x_i-\mu)^2}{2\sigma^2})\times\frac{1}{\sqrt{2\pi}}\exp(-\frac{(\mu-100)^2}{2})\times4(\sigma^2)^{-4}\exp(-\frac{2}{\sigma^2}) \]

となる。

これは\(\sigma, \mu\)について、二変数の確率分布となっているが、 \(\sigma\)に適当な値を入れれば正規分布、\(\mu\)に適当な値を入れれば逆ガンマ分布になる。

この事後分布を用いての予測をすると、

\[ P(X\leq a)=\int_{-\infty}^adx\int \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})p(\mu,\sigma^2\mid D)d\mu d\sigma \]

と計算できる。

自然な共役分布

一般的にはパラメータの事後分布の密度関数を正確に計算したり、 それを用いてパラメータ推定値の平均や分散、あるいは予測分布を計算するのは困難である。

これを克服するための一つの方法は、事前分布として自然共役分布と呼ばれる特定の分布を用いることである。 ベイズ推論の枠組みで、尤度関数はモデリングに依存して決定するものであるが、いくつかの場合にはこの尤度関数と相性のよい事前分布が存在する。

例えば

• 二項分布のパラメータ推定にベータ分布
• 正規分布のパラメータ推定に正規分布、逆ガンマ分布
• ポアソン分布のパタメータ推定にガンマ分布

などがある。 これらを使えば、正確に事後分布を決定することができる。

階層モデル

ある学校で問題数10問のテストを20人の生徒に行ったところ、以下のような結果をえた。 生徒が問題に正解する能力を正解率\(p\)として、 これをパラメータとする二項分布モデルで考えてみよう。

すると、二項分布で想定される分布の様子、あるいは分散の大きさと異なってしまう。

このような場合、正解率に個人差があると考える。

ただし個々の正解率をパラメータ\(p_1,p_2,\ldots,p_{20}\)とするのではなく、 正解率が別の分布に従って決まっていると考える。

このようなモデルを階層モデルという。

マルコフ連鎖モンテカルロ法(MCMC)

例えば階層モデルなどより複雑なモデルなどの場合にも、 コンピュータの助けを借りることでパラメータの推定やデータの予測が可能になる。 これを実現するための一つの方法がMCMCを用いたランダムサンプリングである。

ただし、 この場合には正確な分布を決定するわけではなく、 事後分布に従ったランダムサンプリングが近似的にできるのみであることに注意しよう。

問題設定によっては収束しにくい分布などもあるので、上の自然な共役分布を使った場合と違ってより慎重な運用が必要になる。

MCMCの概要

確率密度関数が与えられたとき、その密度関数に従ってサンプリングを行いたい。 よく知られている確率密度関数であれば、以下のようにRでも簡単にサンプリングすることができる。

N<-10000
x<-rnorm(N)
hist(x)

しかし、より複雑なモデリングを考える場合には、あらかじめ用意されている関数では乱数をサンプリングできない。

モンテカルロ法の例

0または1の値をとり、1をとる確率が円の面積の1/4になるような乱数を生成しよう。 一様分布に従う乱数はあらかじめ用意されているとする。

sample <- function(n){
  x <- c()
  for(i in 1:n){
    r <- runif(2)
    x[i] <- ifelse(r[1]^2+r[2]^2<1, 1, 0)
  }
  return(x)
}

N <- 1000
4*sum(sample(N))/N

## [1] 3.132

メトロポリス法

一様乱数があらかじめ与えられている時、確率密度関数\(f(x)\)に従う乱数\(\{x_1,x_2,\ldots\}\)を以下の手順でサンプリングする。

\(x_i\)から\(x_{i+1}\)を作る方法。

• まず適当な一様乱数から\(\epsilon\)を生成し、次ステップの候補を\(x'=x_i+\epsilon\)とする。
• 次のような条件で実際に\(x_{i+1}=x'\)とするか、あるいは\(x_{i+1}=x_i\)のままにするかを決める。
• \(r=f(x_{i+1})/f(x_i)\)とする。
• \(r\)が1より大きい時、\(x_{i+1}=x'\)とする。
• \(r\)が1より小さい場合、確率\(r\)で\(x_{i+1}=x'\)にする。
• 実際には、0から1の間の一様乱数\(q\)を発生させ、\(r>q\)なら\(x_{i+1}=x'\)とし\(r<q\)なら\(x_{i+1}=x_i\)とする。

p <- function(x){
  return(-6*x*(x-1))
}
x<-0:100/100
plot(x,p(x),'l')

step <- function(x){
  x_next <- x + runif(1, min=-1, max=1)#この幅は調整する
  r <- p(x_next)/p(x)
  if(r > 1){
    return(x_next)
  }
  else{
    q <- runif(1)
    if(r>q) return(x_next)
    else return(x)
  }
}
sample <- c()
sample[1] <- runif(1)
N <- 1000
for(i in 1:(N-1)){
  sample[i+1] <- step(sample[i])
}
hist(sample)

Stanを用いたMCMCサンプリング

Stanの使用法

.stanファイルに分析するdata, parameter, modelを記述する。

実際に分析するデータや、分析結果の処理はRで行う。

実際にstanのコードとそれを用いた推定について、簡単なモデルから順番に見ていこう。

二項分布のパラメータ推定

コインを50回投げて28回表が出たとする。 このデータを元に、コインの表の出る確率\(p\)を推定する。

まずはstanファイルを用意して、以下のように記述する。

data{
  int D;  //成功した回数
  int N;  //試行回数
}
  
parameters{
  real<lower=0,upper=1> p;  //成功確率
}

model{
  D ~ binomial(N,p);  //二項分布モデル
}

次にRのスクリプトに以下のように記述する。

library(rstan)

## Warning: package 'rstan' was built under R version 3.4.3

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.4.4

## Loading required package: StanHeaders

## Warning: package 'StanHeaders' was built under R version 3.4.3

## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

d <- list(D=28, N=50)
fit1 <- stan('binom1.stan', data=d)

## 
## SAMPLING FOR MODEL 'd06ec3f2e0cf24bb2a3cdc911e8284d0' NOW (CHAIN 1).
## 
## Gradient evaluation took 1.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.016615 seconds (Warm-up)
##                0.012988 seconds (Sampling)
##                0.029603 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd06ec3f2e0cf24bb2a3cdc911e8284d0' NOW (CHAIN 2).
## 
## Gradient evaluation took 8e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.016539 seconds (Warm-up)
##                0.017268 seconds (Sampling)
##                0.033807 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd06ec3f2e0cf24bb2a3cdc911e8284d0' NOW (CHAIN 3).
## 
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.015793 seconds (Warm-up)
##                0.012782 seconds (Sampling)
##                0.028575 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd06ec3f2e0cf24bb2a3cdc911e8284d0' NOW (CHAIN 4).
## 
## Gradient evaluation took 4e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.013623 seconds (Warm-up)
##                0.013948 seconds (Sampling)
##                0.027571 seconds (Total)

fit1

## Inference for Stan model: d06ec3f2e0cf24bb2a3cdc911e8284d0.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## p      0.56    0.00 0.07   0.42   0.51   0.56   0.60   0.68  1677    1
## lp__ -36.20    0.02 0.72 -38.16 -36.38 -35.93 -35.75 -35.70  1901    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:35:42 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

推定結果は以下のようにみることができる。

パラメータのサンプリング、横軸がステップ数

stan_trace(fit1,pars="p")

パラメータのサンプルのヒストグラム

stan_hist(fit1,pars="p")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit1,pars="p",separate_chains = TRUE)

パラメータの自己相関

stan_ac(fit1,pars = "p",separate_chains = TRUE)

収束の判断。

• Rhatが1に近い（1.05以下が一つの目安）
• チェインごとの確率密度が重なっている
• 自己相関が低い。定常分布に収束したら相関はないはず

サンプリングした上で平均など色々な関数の値を計算することができる。 またそれを用いて未知のデータについて予測できる。

例えば次にコインを投げて表の出る確率は、以下のように予測できる。

p <- rstan::extract(fit1)$p
mean(p)

## [1] 0.5565015

ポアソン分布のパラメータ推定

次にポアソン分布のパラメータ推定をしてみよう。

前に扱った事故件数のデータを用いる。

まず.stanファイルにモデルについて記述する。

data{
  int N;  //サンプル数
  int x[N];  //観測値
}
  
parameters{
  real<lower=0> lambda;  //ポアソン分布のパラメータ
}

model{
  for (i in 1:N){
    x[i] ~ poisson(lambda);  //ポアソン分布モデル
    }
}

次にデータを用意し、Rで分析する。

x <- c(8,4,5,5,7,9,7,5,5,4)
d <- list(x=x, N=length(x))
fit2 <- stan('poisson1.stan', data=d)

## 
## SAMPLING FOR MODEL 'd9d5669abb02a1e798ca9ce33cd34eb4' NOW (CHAIN 1).
## 
## Gradient evaluation took 1.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.024758 seconds (Warm-up)
##                0.02992 seconds (Sampling)
##                0.054678 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd9d5669abb02a1e798ca9ce33cd34eb4' NOW (CHAIN 2).
## 
## Gradient evaluation took 1e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.018442 seconds (Warm-up)
##                0.017348 seconds (Sampling)
##                0.03579 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd9d5669abb02a1e798ca9ce33cd34eb4' NOW (CHAIN 3).
## 
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.016558 seconds (Warm-up)
##                0.013974 seconds (Sampling)
##                0.030532 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'd9d5669abb02a1e798ca9ce33cd34eb4' NOW (CHAIN 4).
## 
## Gradient evaluation took 5e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.014811 seconds (Warm-up)
##                0.013201 seconds (Sampling)
##                0.028012 seconds (Total)

fit2

## Inference for Stan model: d9d5669abb02a1e798ca9ce33cd34eb4.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## lambda  6.02    0.02 0.78  4.61  5.46  5.98  6.52  7.71  1418    1
## lp__   47.00    0.02 0.71 45.00 46.83 47.27 47.45 47.50  1517    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:36:44 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

パラメータのサンプリング、横軸がステップ数

stan_trace(fit2,pars="lambda")

パラメータのサンプルのヒストグラム

stan_hist(fit2,pars="lambda")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit2,pars="lambda",separate_chains = TRUE)

パラメータの自己相関

stan_ac(fit2,pars = "lambda",separate_chains = TRUE)

事故件数の平均の予測値は以下のようになる。

lambda <- rstan::extract(fit2)$lambda
mean(lambda)

## [1] 6.016857

正規分布パラメータの推定

次に正規分布のパラメータを推定してみよう。

前と同じ、製品の規格検査のデータを使おう。

data{
  int N;      //サンプル数
  real y[N];  //観測値
}

parameters{
  real mu;    //平均
  real<lower=0> sigma; //標準偏差
}

model{
  mu ~ normal(0,100);   //平均の事前分布
  sigma ~ cauchy(0,5);  //標準偏差の事前分布
  for(i in 1:N){
    y[i] ~ normal(mu,sigma);
  }
}

y <- c(102.3,100.4,99.3,100.2,100.3,99.4,101.5,99.3,98.9,100.1)
d <- list(y=y,N=length(y))
fit3 <- stan(file = "normal1.stan", data = d)

## 
## SAMPLING FOR MODEL 'c81f4001bf9438379267dee6a2a4f8a2' NOW (CHAIN 1).
## 
## Gradient evaluation took 1.5e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.026445 seconds (Warm-up)
##                0.018703 seconds (Sampling)
##                0.045148 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'c81f4001bf9438379267dee6a2a4f8a2' NOW (CHAIN 2).
## 
## Gradient evaluation took 5e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.022583 seconds (Warm-up)
##                0.018683 seconds (Sampling)
##                0.041266 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'c81f4001bf9438379267dee6a2a4f8a2' NOW (CHAIN 3).
## 
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.020552 seconds (Warm-up)
##                0.030861 seconds (Sampling)
##                0.051413 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'c81f4001bf9438379267dee6a2a4f8a2' NOW (CHAIN 4).
## 
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.022397 seconds (Warm-up)
##                0.028378 seconds (Sampling)
##                0.050775 seconds (Total)

fit3

## Inference for Stan model: c81f4001bf9438379267dee6a2a4f8a2.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean   sd  2.5%   25%    50%    75%  97.5% n_eff Rhat
## mu    100.16    0.01 0.40 99.36 99.92 100.16 100.40 100.92  1899    1
## sigma   1.23    0.01 0.34  0.75  0.99   1.16   1.39   2.03  1891    1
## lp__   -6.63    0.03 1.08 -9.44 -7.02  -6.31  -5.87  -5.60  1249    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:37:36 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

パラメータのサンプリング、横軸がステップ数

stan_trace(fit3, pars=c("mu","sigma"))

パラメータのサンプルのヒストグラム

stan_hist(fit3,pars=c("mu","sigma"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit3, pars=c("mu","sigma"),separate_chains = TRUE)

パラメータの自己相関

stan_ac(fit3,pars = c("mu","sigma"),separate_chains = TRUE)

t検定

上の正規分布のパラメータの推定の応用で、 二つのグループの母平均を比較してみよう。

次のように、二つのクラスでテストをした結果が得られたとする。 このニクラスの平均点に差があると言えるかどうか、ベイズ推定しよう。

各クラスの点数は正規分布に従うと仮定し、 そのパラメータをそれぞれmu_x, sigma_xとmu_y, sigma_yとしてモデルを作る。

data{
  int N;  //サンプル数
  real y[N];  //群1
  real x[N];  //群2
}

parameters{
  real mu_y;  //yの平均値
  real<lower=0> sigma_y;  //yの標準偏差
  real mu_x;  //xの平均値
  real<lower=0> sigma_x;  //xの標準偏差
}

model{
  mu_y ~ normal(0,100); //無情報事前分布
  sigma_y ~ cauchy(0,5);//無情報事前分布
  mu_x ~ normal(0,100); //無情報事前分布
  sigma_x ~ cauchy(0,5);//無情報事前分布
  
  y ~ normal(mu_y,sigma_y);
  x ~ normal(mu_x,sigma_x);
}

generated quantities{
  real diff;  //2群の平均の差
  diff = mu_x-mu_y;
}

library(rstan)
N<-10
y<-c(60,39,75,64,36,45,38,28,3,9)
x<-c(73,67,19,100,100,0,58,39,84,62)

d<-list(N=N, x=x, y=y)
fit <- stan(file = "ttest.stan", data = d)

## 
## SAMPLING FOR MODEL 'b6f16321d07933d8582b0bc79e69a9e1' NOW (CHAIN 1).
## 
## Gradient evaluation took 1.5e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.05534 seconds (Warm-up)
##                0.023017 seconds (Sampling)
##                0.078357 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'b6f16321d07933d8582b0bc79e69a9e1' NOW (CHAIN 2).
## 
## Gradient evaluation took 7e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.049284 seconds (Warm-up)
##                0.021569 seconds (Sampling)
##                0.070853 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'b6f16321d07933d8582b0bc79e69a9e1' NOW (CHAIN 3).
## 
## Gradient evaluation took 5e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.04901 seconds (Warm-up)
##                0.032627 seconds (Sampling)
##                0.081637 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'b6f16321d07933d8582b0bc79e69a9e1' NOW (CHAIN 4).
## 
## Gradient evaluation took 6e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.047897 seconds (Warm-up)
##                0.022158 seconds (Sampling)
##                0.070055 seconds (Total)

fit

## Inference for Stan model: b6f16321d07933d8582b0bc79e69a9e1.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##           mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu_y     39.42    0.13  7.58  24.04  34.69  39.41  44.25  54.70  3244    1
## sigma_y  23.56    0.11  6.02  14.88  19.44  22.51  26.45  38.07  2890    1
## mu_x     59.48    0.20 11.02  37.22  52.83  59.43  66.14  81.59  3173    1
## sigma_x  34.13    0.18  8.91  21.72  27.87  32.56  38.53  55.64  2400    1
## diff     20.06    0.24 13.57  -6.70  11.51  19.90  28.86  47.54  3207    1
## lp__    -77.80    0.04  1.56 -81.70 -78.60 -77.49 -76.63 -75.83  1548    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:38:36 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

diff <- rstan::extract(fit)$diff
p <- sum(ifelse(diff>0,1,0))/length(diff)
p

## [1] 0.939

plot(density(diff))

パラメータのサンプリング、横軸がステップ数

stan_trace(fit)

パラメータのサンプルのヒストグラム

stan_hist(fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit,separate_chains = TRUE)

パラメータの自己相関

stan_ac(fit3,separate_chains = TRUE)

線形回帰

回帰分析もベイズ推定の枠組みで扱うことができる。 ここでは以下のデータを用いて単回帰分析を行う。

data{
  int N;  //サンプル数
  real x[N];  //説明変数
  real y[N];  //目的変数
}

parameters{
  real a; //傾き
  real b; //切片
  real<lower=0> sigma;  //誤差の標準偏差
}

model{
  for (i in 1:N){
    y[i] ~ normal(b+a*x[i],sigma);
  }
}

library(rstan)
N<-30
x<-runif(N)
noise<-rnorm(N,mean=0,sd=0.1)
a<-1
b<-2
y<-b+a*x+noise
d<-list(x=x, y=y, N=N)
fit<-stan(file='linreg.stan',data=d)

## 
## SAMPLING FOR MODEL '588538e1cebbf0d9aad8232558d90dcc' NOW (CHAIN 1).
## 
## Gradient evaluation took 1.7e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.056297 seconds (Warm-up)
##                0.055837 seconds (Sampling)
##                0.112134 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '588538e1cebbf0d9aad8232558d90dcc' NOW (CHAIN 2).
## 
## Gradient evaluation took 8e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.05539 seconds (Warm-up)
##                0.058719 seconds (Sampling)
##                0.114109 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '588538e1cebbf0d9aad8232558d90dcc' NOW (CHAIN 3).
## 
## Gradient evaluation took 1.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.056301 seconds (Warm-up)
##                0.053084 seconds (Sampling)
##                0.109385 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '588538e1cebbf0d9aad8232558d90dcc' NOW (CHAIN 4).
## 
## Gradient evaluation took 1.1e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.062233 seconds (Warm-up)
##                0.059078 seconds (Sampling)
##                0.121311 seconds (Total)

fit

## Inference for Stan model: 588538e1cebbf0d9aad8232558d90dcc.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## a      0.88    0.00 0.08  0.72  0.82  0.88  0.93  1.02  1963    1
## b      2.06    0.00 0.04  1.99  2.04  2.06  2.09  2.13  1914    1
## sigma  0.11    0.00 0.02  0.09  0.10  0.11  0.12  0.15  2253    1
## lp__  48.59    0.03 1.24 45.51 48.01 48.86 49.50 50.05  1767    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:39:26 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

パラメータのサンプリング、横軸がステップ数

stan_trace(fit)

パラメータのサンプルのヒストグラム

stan_hist(fit)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit,separate_chains = TRUE)

パラメータの自己相関

stan_ac(fit3,separate_chains = TRUE)

階層ベイズモデル

統計的モデリングの概要

データを発生させるメカニズムがある確率分布で記述できるとし、 それをよく説明する確率分布を作る。

階層ベイズとは

データが単純に一つの確率分布から発生するのではなく、 パラメータが発生するメカニズムとしてさらに確率分布を用いるなど、 いくつかの確率分布を積み上げてモデルを作る。

例えば個体差、場所差などを入れることができる。

具体的に、以下の事例を見ながら、階層モデリングの様子を把握しよう。

ロジスティック回帰

ある20人の生徒に小テストを行った。 テストは10問からなり、結果は以下のようになった。

x <- c(10,9,6,7,8,3,4,9,3,7,10,2,0,5,6,3,1,9,0,0)

まずは二項分布モデルに基づいて、生徒の学力を正解率\(p\)として推定しよう。 尤度関数を計算すると、以下のようになる。 \[ L(p) = C~p^{102}(1-p)^{98} \] ここで\(C\)は適当な定数。

事前分布を一様分布として、事後分布は上で求めた結果からベータ分布\(B(103,99)\)となる。 このとき\(p\)の平均はおよそ0.51で、これに対する二項分布の分散はおよそ2.5である。

一方でデータから計算できる分散は12.1となる。

このように二項分布モデルで推定した分散より、 実際のデータの分散が大きくなる現象を過分散という。

これを解消するために、階層モデルを用いて予測する。

例えばここでは、生徒の個人差のパラメータ\(r_i\)を用いることにしよう。 またこの個人差のパラメータは標準偏差sigmaの正規分布の従うとする。

data{
  int N;  //サンプル数
  int M;  //テスト問題数
  int x[N]; //正解数
}

parameters{
  real r[N];  //個人差
  real beta;  //共通パラメータ
}

transformed parameters{
  real<lower=0,upper=1> p[N]; //個人の正解率
  for (i in 1:N){
    p[i] = inv_logit(beta+r[i]);
  }
}

model{
  for (i in 1:N){
    x[i] ~ binomial(M,p[i]);
  }
}

上のコードを実行する。

library(rstan)
N<-20
x <- c(10,9,6,7,8,3,4,9,3,7,10,2,0,5,6,3,1,9,0,0)
d<-list(x=x, N=N, M=10)
fit<-stan(file='logreg.stan',data=d)

## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config/compiler/clang.hpp:200:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##           ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
##         ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Core:531:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/operator_unary_plus.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/scal/fun/constants.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/constants/constants.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/tools/convert_from_string.hpp:15:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast.hpp:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast/try_lexical_convert.hpp:42:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast/detail/converter_lexical.hpp:52:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/container/container_fwd.hpp:61:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/container/detail/std_fwd.hpp:27:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/__config:390:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d76993776b.cpp:487:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rstan/include/rstan/rstaninc.hpp:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rstan/include/rstan/stan_fit.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/services/optimize/bfgs.hpp:11:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/optimization/bfgs.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/optimization/lbfgs_update.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/circular_buffer.hpp:54:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/circular_buffer/details.hpp:20:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/move.hpp:30:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/iterator.hpp:27:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/iterator_traits.hpp:29:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/__config:390:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
##     static void set_zero_all_adjoints() {
##                 ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:45:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
##     static void set_zero_all_adjoints_nested() {
##                 ^
## In file included from filea3d76993776b.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
##       size_t fft_next_good_size(size_t N) {
##              ^
## 19 warnings generated.
## 
## SAMPLING FOR MODEL 'logreg' NOW (CHAIN 1).
## 
## Gradient evaluation took 2.5e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.144554 seconds (Warm-up)
##                0.135262 seconds (Sampling)
##                0.279816 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'logreg' NOW (CHAIN 2).
## 
## Gradient evaluation took 1.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.158835 seconds (Warm-up)
##                0.138877 seconds (Sampling)
##                0.297712 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'logreg' NOW (CHAIN 3).
## 
## Gradient evaluation took 1.3e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.140096 seconds (Warm-up)
##                0.132215 seconds (Sampling)
##                0.272311 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'logreg' NOW (CHAIN 4).
## 
## Gradient evaluation took 1.1e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.147344 seconds (Warm-up)
##                0.125105 seconds (Sampling)
##                0.272449 seconds (Total)

fit

## Inference for Stan model: logreg.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff
## sigma    2.45    0.02 0.67    1.44    1.99    2.35    2.80    4.10  1851
## r[1]     3.40    0.03 1.58    0.90    2.30    3.20    4.26    7.01  2265
## r[2]     2.10    0.03 1.12    0.14    1.31    2.04    2.79    4.44  1889
## r[3]     0.41    0.02 0.86   -1.30   -0.15    0.42    0.97    2.11  1384
## r[4]     0.85    0.02 0.89   -0.81    0.26    0.84    1.44    2.63  1448
## r[5]     1.39    0.02 0.95   -0.35    0.74    1.34    1.99    3.37  1718
## r[6]    -0.86    0.02 0.93   -2.83   -1.45   -0.83   -0.24    0.89  1622
## r[7]    -0.42    0.02 0.88   -2.23   -0.97   -0.42    0.17    1.29  1544
## r[8]     2.14    0.03 1.16    0.14    1.37    2.02    2.81    4.79  1982
## r[9]    -0.85    0.02 0.89   -2.67   -1.41   -0.85   -0.26    0.87  1464
## r[10]    0.85    0.03 0.90   -0.83    0.24    0.84    1.45    2.68  1295
## r[11]    3.41    0.03 1.61    0.88    2.29    3.16    4.27    7.25  2234
## r[12]   -1.37    0.02 0.96   -3.36   -1.98   -1.34   -0.73    0.40  1598
## r[13]   -3.37    0.03 1.53   -6.87   -4.20   -3.15   -2.32   -0.96  2593
## r[14]    0.00    0.02 0.86   -1.65   -0.56   -0.01    0.58    1.70  1295
## r[15]    0.41    0.02 0.87   -1.31   -0.15    0.41    0.98    2.14  1419
## r[16]   -0.85    0.02 0.91   -2.62   -1.45   -0.85   -0.25    0.88  1434
## r[17]   -2.11    0.02 1.11   -4.64   -2.72   -2.02   -1.36   -0.16  1993
## r[18]    2.10    0.02 1.11    0.08    1.33    2.04    2.76    4.52  2049
## r[19]   -3.44    0.03 1.63   -7.29   -4.30   -3.22   -2.32   -0.91  2362
## r[20]   -3.42    0.03 1.54   -7.10   -4.26   -3.22   -2.36   -0.91  2367
## beta     0.01    0.02 0.63   -1.25   -0.38    0.01    0.38    1.24   817
## p[1]     0.94    0.00 0.07    0.75    0.91    0.96    0.99    1.00  4000
## p[2]     0.86    0.00 0.10    0.62    0.80    0.88    0.94    0.99  4000
## p[3]     0.59    0.00 0.14    0.31    0.50    0.60    0.69    0.84  4000
## p[4]     0.69    0.00 0.13    0.41    0.60    0.69    0.79    0.90  4000
## p[5]     0.78    0.00 0.12    0.51    0.70    0.79    0.87    0.96  4000
## p[6]     0.32    0.00 0.14    0.08    0.21    0.30    0.41    0.63  4000
## p[7]     0.41    0.00 0.15    0.15    0.30    0.40    0.51    0.70  4000
## p[8]     0.86    0.00 0.10    0.62    0.81    0.88    0.94    0.99  4000
## p[9]     0.32    0.00 0.13    0.10    0.22    0.30    0.40    0.61  4000
## p[10]    0.68    0.00 0.14    0.39    0.59    0.70    0.78    0.91  4000
## p[11]    0.94    0.00 0.07    0.75    0.91    0.96    0.99    1.00  4000
## p[12]    0.23    0.00 0.12    0.05    0.14    0.21    0.30    0.50  4000
## p[13]    0.06    0.00 0.07    0.00    0.02    0.04    0.09    0.24  4000
## p[14]    0.50    0.00 0.14    0.23    0.40    0.50    0.60    0.78  4000
## p[15]    0.59    0.00 0.14    0.30    0.49    0.60    0.69    0.85  4000
## p[16]    0.32    0.00 0.14    0.09    0.21    0.31    0.41    0.61  4000
## p[17]    0.14    0.00 0.10    0.01    0.07    0.12    0.19    0.38  4000
## p[18]    0.86    0.00 0.10    0.62    0.81    0.88    0.93    0.99  4000
## p[19]    0.06    0.00 0.07    0.00    0.01    0.04    0.09    0.26  4000
## p[20]    0.06    0.00 0.07    0.00    0.01    0.04    0.09    0.25  4000
## lp__  -116.98    0.11 4.00 -125.92 -119.39 -116.61 -114.10 -110.30  1276
##       Rhat
## sigma 1.00
## r[1]  1.00
## r[2]  1.00
## r[3]  1.00
## r[4]  1.00
## r[5]  1.00
## r[6]  1.00
## r[7]  1.00
## r[8]  1.00
## r[9]  1.00
## r[10] 1.00
## r[11] 1.00
## r[12] 1.00
## r[13] 1.00
## r[14] 1.00
## r[15] 1.00
## r[16] 1.00
## r[17] 1.00
## r[18] 1.00
## r[19] 1.00
## r[20] 1.00
## beta  1.01
## p[1]  1.00
## p[2]  1.00
## p[3]  1.00
## p[4]  1.00
## p[5]  1.00
## p[6]  1.00
## p[7]  1.00
## p[8]  1.00
## p[9]  1.00
## p[10] 1.00
## p[11] 1.00
## p[12] 1.00
## p[13] 1.00
## p[14] 1.00
## p[15] 1.00
## p[16] 1.00
## p[17] 1.00
## p[18] 1.00
## p[19] 1.00
## p[20] 1.00
## lp__  1.00
## 
## Samples were drawn using NUTS(diag_e) at Sat Apr 20 06:40:58 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

パラメータのサンプリング、横軸がステップ数

stan_trace(fit)

## 'pars' not specified. Showing first 10 parameters by default.

パラメータのサンプルのヒストグラム

stan_hist(fit)

## 'pars' not specified. Showing first 10 parameters by default.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

パラメータ事後分布の密度推定

stan_dens(fit)

## 'pars' not specified. Showing first 10 parameters by default.

パラメータの自己相関

stan_ac(fit)

## 'pars' not specified. Showing first 10 parameters by default.

各生徒の正解率の推定は以下のようになる。

stan_hist(fit, pars = "p")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

状態空間モデル

状態方程式と観測方程式

観測値\(y_n\)を直接説明するのではなく、真の状態\(x_n\)の変化とその観測に分けて考える。

真の状態の様子を説明する式をシステムモデルや状態方程式と呼ぶ。 \[ x_n = F_n(x_{n-1}) + G_n(v_n) \] ここで\(v_n\)はシステムノイズといわれ、ある確率分布に従う乱数。

また真の状態から観測値が発生する様子を説明する式を観測モデルや観測方程式という。 \[ y_n = H_n(x_n) + w_n \] ここで\(w_n\)は観測ノイズといわれ、ある確率分布に従う乱数。

一番基本的な例は線形ガウス型モデルと言われるもので、 上の式において、線形関数\(F_n(x)=F_nx, G_n(v)=G_nv, H_n(x)=H_nx\)とし、\(v_n\sim N(0,Q_n), w_n\sim N(0,R_n)\)は正規分布に従う乱数。

また時系列モデルの基本的な例であるARモデル

\[ y_n = a_1y_{n-1}+a_2y_{n-2}+v_n,~v_n\sim N(0,\sigma^2) \]

は、状態空間モデルの言葉を用いると、 \[ x_n=\begin{pmatrix}y_n\\y_{n-1}\end{pmatrix}, F=\begin{pmatrix}a_1&a_2\\1&0\end{pmatrix}, G=\begin{pmatrix}1\\0\end{pmatrix},H=\begin{pmatrix}1&0\end{pmatrix}, Q=\sigma^2, R=0 \]

として表すことができる。

ローカルレベルモデル

ローカルレベルモデルは次のような状態方程式と観測方程式をもつ。

tを時刻、\(x_t\)を状態、\(y_t\)を観測値として、 \[ y_t=x_t+\epsilon \] \[ x_t=x_{t-1}+\delta \] \[ \epsilon \sim N(0,\sigma_1) \] \[ \delta \sim N(0,\sigma_2) \]

状態方程式はいわゆるランダムウォークになっている。

z<-rnorm(100)
y<-0
for (i in 1:100){
  y[i+1]<-y[i]+z[i]
}
plot(y,type='l')

推定するパラメータは状態\(x_t\)及び誤差の分散 \(\sigma_1,\sigma_2\)

data {
  int T; //サンプルサイズ
  real y[T]; //観測値
}

parameters {
  real x[T]; //状態
  real<lower=0> sigma_s; //状態のノイズの分散
  real<lower=0> sigma_o; //観測のノイズの分散
}

model {
  //観測方程式
  for (i in 1:T){
    y[i] ~ normal(x[i], sigma_o);
  }
  //状態方程式
  for (i in 2:T){
    x[i] ~ normal(x[i-1], sigma_s);
  }
}

#ナイル川データ
Nile

## Time Series:
## Start = 1871 
## End = 1970 
## Frequency = 1 
##   [1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935 1110  994
##  [15] 1020  960 1180  799  958 1140 1100 1210 1150 1250 1260 1220 1030 1100
##  [29]  774  840  874  694  940  833  701  916  692 1020 1050  969  831  726
##  [43]  456  824  702 1120 1100  832  764  821  768  845  864  862  698  845
##  [57]  744  796 1040  759  781  865  845  944  984  897  822 1010  771  676
##  [71]  649  846  812  742  801 1040  860  874  848  890  744  749  838 1050
##  [85]  918  986  797  923  975  815 1020  906  901 1170  912  746  919  718
##  [99]  714  740

y<-list(y=as.numeric(Nile),T=length(Nile))
y

## $y
##   [1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935 1110  994
##  [15] 1020  960 1180  799  958 1140 1100 1210 1150 1250 1260 1220 1030 1100
##  [29]  774  840  874  694  940  833  701  916  692 1020 1050  969  831  726
##  [43]  456  824  702 1120 1100  832  764  821  768  845  864  862  698  845
##  [57]  744  796 1040  759  781  865  845  944  984  897  822 1010  771  676
##  [71]  649  846  812  742  801 1040  860  874  848  890  744  749  838 1050
##  [85]  918  986  797  923  975  815 1020  906  901 1170  912  746  919  718
##  [99]  714  740
## 
## $T
## [1] 100

library(rstan)
fit <- stan(file = 'locallevel.stan', data = y,
                 iter = 4000, chains = 4)

## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config.hpp:39:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/config/compiler/clang.hpp:200:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##           ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
##         ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Core:531:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/operator_unary_plus.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/scal/fun/constants.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/constants/constants.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/math/tools/convert_from_string.hpp:15:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast.hpp:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast/try_lexical_convert.hpp:42:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/lexical_cast/detail/converter_lexical.hpp:52:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/container/container_fwd.hpp:61:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/container/detail/std_fwd.hpp:27:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/__config:390:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filea3d72412919.cpp:423:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rstan/include/rstan/rstaninc.hpp:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rstan/include/rstan/stan_fit.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/services/optimize/bfgs.hpp:11:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/optimization/bfgs.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/optimization/lbfgs_update.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/circular_buffer.hpp:54:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/circular_buffer/details.hpp:20:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/move.hpp:30:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/iterator.hpp:27:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/iterator_traits.hpp:29:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/__config:390:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
##     static void set_zero_all_adjoints() {
##                 ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:45:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints_nested.hpp:17:17: warning: 'static' function 'set_zero_all_adjoints_nested' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
##     static void set_zero_all_adjoints_nested() {
##                 ^
## In file included from filea3d72412919.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:58:
## /Library/Frameworks/R.framework/Versions/3.4/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:17:14: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
##       size_t fft_next_good_size(size_t N) {
##              ^
## 19 warnings generated.
## 
## SAMPLING FOR MODEL 'locallevel' NOW (CHAIN 1).
## 
## Gradient evaluation took 3.9e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 2.51381 seconds (Warm-up)
##                1.55357 seconds (Sampling)
##                4.06738 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'locallevel' NOW (CHAIN 2).
## 
## Gradient evaluation took 1.9e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 5.50166 seconds (Warm-up)
##                1.43336 seconds (Sampling)
##                6.93502 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'locallevel' NOW (CHAIN 3).
## 
## Gradient evaluation took 1.6e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 3.06602 seconds (Warm-up)
##                1.43698 seconds (Sampling)
##                4.50301 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'locallevel' NOW (CHAIN 4).
## 
## Gradient evaluation took 2.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 2.56642 seconds (Warm-up)
##                1.81625 seconds (Sampling)
##                4.38267 seconds (Total)

## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

#状態の事後分布
x<-rstan::extract(fit)$x
#状態の推定結果を可視化する
head(x)

##           
## iterations     [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
##       [1,] 1236.779 1173.922 1216.205 1173.419 1205.783 1094.579 1117.453
##       [2,] 1061.119 1073.747 1073.192 1092.841 1120.160 1151.270 1160.773
##       [3,] 1131.055 1081.093 1067.785 1073.910 1117.070 1094.862 1059.443
##       [4,] 1170.974 1193.193 1182.012 1115.981 1161.813 1189.283 1185.489
##       [5,] 1181.831 1134.036 1076.515 1051.846 1122.995 1155.547 1155.634
##       [6,] 1061.359 1046.477 1031.731 1037.765 1016.486 1018.686 1028.022
##           
## iterations     [,8]     [,9]    [,10]     [,11]    [,12]    [,13]    [,14]
##       [1,] 1088.280 1196.890 1191.363 1123.8066 1039.762 1082.412 1055.406
##       [2,] 1169.424 1167.613 1106.015 1088.0244 1050.207 1071.467 1062.692
##       [3,] 1023.126 1024.421 1004.681  998.6164 1059.085 1045.882 1006.567
##       [4,] 1161.110 1159.429 1140.953 1103.5148 1066.596 1064.717 1061.555
##       [5,] 1173.486 1178.879 1057.277 1102.2093 1034.894 1068.786 1044.080
##       [6,] 1048.336 1077.823 1053.784 1015.7844 1002.840 1017.314 1013.813
##           
## iterations    [,15]     [,16]    [,17]    [,18]    [,19]    [,20]    [,21]
##       [1,] 1097.750 1119.6450 1118.429 1159.762 1154.654 1113.724 1047.321
##       [2,] 1064.359 1050.2200 1054.509 1024.183 1009.167 1096.831 1120.052
##       [3,] 1041.191 1069.9576 1097.492 1034.108 1057.742 1071.913 1097.509
##       [4,] 1062.316 1047.3673 1098.449 1076.825 1086.436 1058.554 1052.262
##       [5,] 1068.409 1035.1385 1077.090 1049.836 1057.128 1091.022 1020.879
##       [6,] 1007.994  981.1574 1034.272 1033.035 1033.235 1040.687 1048.750
##           
## iterations    [,22]    [,23]    [,24]    [,25]    [,26]    [,27]     [,28]
##       [1,] 1121.155 1144.588 1189.702 1164.683 1092.801 1039.885 1024.2595
##       [2,] 1154.453 1097.230 1062.160 1046.596 1145.282 1118.766 1074.9030
##       [3,] 1075.032 1089.237 1073.210 1072.244 1017.636 1003.850  955.9608
##       [4,] 1097.042 1073.283 1135.807 1106.210 1061.959 1006.640  964.3959
##       [5,] 1074.298 1132.300 1137.044 1081.825 1087.886 1044.645  978.7098
##       [6,] 1081.478 1077.775 1109.390 1094.758 1046.645 1005.908  960.4691
##           
## iterations     [,29]    [,30]    [,31]    [,32]    [,33]    [,34]    [,35]
##       [1,]  897.0412 836.1469 797.8580 872.2364 810.7340 847.0837 857.5340
##       [2,] 1011.4522 963.9545 844.7967 787.6629 856.3940 822.7337 875.2733
##       [3,]  918.1668 928.7545 934.7954 857.3411 880.2910 861.3335 847.9220
##       [4,]  935.7012 892.4798 894.9718 925.9180 877.9715 849.5332 815.0384
##       [5,]  954.0886 887.8112 842.8156 762.8177 819.4247 800.5792 807.1217
##       [6,]  928.9708 911.5253 903.2390 914.9737 916.4301 894.2090 884.8042
##           
## iterations    [,36]    [,37]    [,38]    [,39]    [,40]    [,41]    [,42]
##       [1,] 888.0906 799.2115 919.1464 928.3189 833.6906 775.9138 759.0595
##       [2,] 906.0097 843.9290 791.7238 841.1972 863.0524 813.4629 789.0172
##       [3,] 839.7131 888.4171 871.9712 847.9283 833.2919 835.5532 830.1818
##       [4,] 843.3867 875.1907 887.9921 867.4470 862.6988 823.8933 827.4089
##       [5,] 882.7364 908.2602 886.0696 911.0613 916.8791 917.2998 914.9296
##       [6,] 862.4021 813.0055 832.1566 836.9579 824.0952 794.6144 775.6408
##           
## iterations    [,43]    [,44]    [,45]    [,46]    [,47]    [,48]    [,49]
##       [1,] 791.2762 816.6553 902.0062 934.6704 781.9926 844.6451 778.2953
##       [2,] 775.3895 824.9568 919.8737 908.9890 945.9218 929.9072 877.1298
##       [3,] 819.7442 886.3627 891.1624 924.8808 933.7311 883.7029 836.4440
##       [4,] 825.2323 849.6089 875.2244 899.3448 906.1389 919.5251 928.5876
##       [5,] 836.1695 827.3422 743.6148 796.8305 903.1979 853.2184 895.6513
##       [6,] 738.0358 766.0654 810.1914 815.5965 816.0819 794.9130 845.1757
##           
## iterations    [,50]    [,51]    [,52]    [,53]    [,54]    [,55]    [,56]
##       [1,] 823.5067 861.5126 879.2801 882.6511 817.6814 790.8407 882.4869
##       [2,] 902.4336 856.1408 867.1598 880.2548 901.0255 884.0601 882.2272
##       [3,] 808.9682 791.2732 829.2461 894.1141 840.0467 849.6578 877.0869
##       [4,] 876.9843 871.2664 852.3374 820.4307 774.0284 806.4070 772.0719
##       [5,] 813.5792 799.4080 838.9480 878.2779 894.2561 841.2829 767.0461
##       [6,] 855.5766 853.9066 858.8412 884.2577 886.1742 895.6434 834.0894
##           
## iterations    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]    [,63]
##       [1,] 789.3003 827.3911 902.5914 823.0189 793.8166 784.7088 849.6820
##       [2,] 795.7119 818.8597 775.2390 804.8333 845.3313 858.3280 899.1995
##       [3,] 892.3655 902.4222 921.6823 903.4130 895.4445 899.3707 915.9702
##       [4,] 806.2271 845.4123 841.7912 891.4192 912.7294 905.9324 912.2308
##       [5,] 745.5575 849.7579 942.9978 847.1491 811.7491 829.1953 797.3766
##       [6,] 851.2853 851.6795 784.2314 783.7007 809.6692 794.4852 827.3438
##           
## iterations    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
##       [1,] 970.1539 885.0055 833.3629 903.0856 861.2760 807.5961 692.9054
##       [2,] 958.4216 958.5025 871.6393 888.2539 876.4597 868.2688 838.2963
##       [3,] 892.9297 868.3459 823.5790 843.7150 864.8400 851.5318 811.6447
##       [4,] 898.2956 887.0869 837.5130 810.3816 880.6576 852.3471 826.2160
##       [5,] 795.2952 822.9157 825.8074 804.3090 847.7225 807.7867 831.3405
##       [6,] 851.6587 862.3359 859.0870 857.4665 871.0682 851.9613 857.6014
##           
## iterations    [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
##       [1,] 732.2878 854.9998 846.0535 859.5370 839.0419 854.8492 805.3545
##       [2,] 781.6948 672.4091 716.3512 798.6791 814.0263 807.0942 784.3869
##       [3,] 841.4866 767.4765 809.1293 804.4957 817.2010 861.6685 881.6933
##       [4,] 879.6754 886.5599 916.0006 938.3529 933.8732 939.9398 946.0245
##       [5,] 827.3608 820.3219 813.8029 830.8050 799.8541 808.9422 871.3972
##       [6,] 831.2107 857.6444 841.1075 852.4603 874.8621 842.1099 830.2392
##           
## iterations    [,78]    [,79]    [,80]    [,81]    [,82]    [,83]    [,84]
##       [1,] 857.7076 860.3977 866.7235 915.5063 932.3734 905.4221 877.6013
##       [2,] 785.4076 782.2046 834.0320 804.2984 783.1051 843.0018 832.3008
##       [3,] 885.0079 851.1942 853.9880 881.1024 843.4134 832.5886 862.7653
##       [4,] 913.3189 872.5768 853.1818 829.4331 819.5214 902.7051 907.6142
##       [5,] 807.7262 852.0669 870.5896 902.3644 887.5319 927.9446 993.5172
##       [6,] 807.2768 806.9863 798.3613 788.6315 797.5077 783.3244 792.2684
##           
## iterations    [,85]    [,86]    [,87]     [,88]    [,89]    [,90]    [,91]
##       [1,] 925.1845 950.3696 951.6005 1014.5602 931.0379 977.9835 892.7869
##       [2,] 867.4759 866.5068 918.2790  886.4097 875.4924 896.4439 867.4428
##       [3,] 909.5484 879.6619 849.0620  892.9505 907.2649 937.0462 962.7322
##       [4,] 884.9065 920.3843 943.3332  937.9837 939.0068 945.2819 974.5606
##       [5,] 901.4785 908.0050 896.1901  881.3720 857.1367 862.6414 880.6562
##       [6,] 790.6572 800.7291 836.1715  874.6512 904.6334 945.0666 890.5244
##           
## iterations    [,92]     [,93]     [,94]    [,95]    [,96]    [,97]
##       [1,] 870.8145 1017.8928 1002.9777 988.6650 895.7644 908.2953
##       [2,] 909.9672  896.2640  881.9934 802.6634 802.0554 824.6263
##       [3,] 985.0146 1017.1438 1012.9133 958.8140 930.9640 898.6679
##       [4,] 934.9321  933.7034  915.9060 944.0167 885.2999 877.9207
##       [5,] 920.9861  917.4130  940.2740 914.5544 843.5148 792.8937
##       [6,] 902.9390  925.8799  973.0033 933.6885 910.5727 865.9566
##           
## iterations    [,98]    [,99]   [,100]
##       [1,] 830.4913 769.7880 777.7221
##       [2,] 755.1011 702.5823 648.3773
##       [3,] 850.6837 857.9988 837.7369
##       [4,] 883.3743 853.0107 884.0148
##       [5,] 835.1081 811.0035 823.4085
##       [6,] 805.2049 816.5582 806.6083

d2<-data.frame(x)
up<-c()
lo<-c()
M<-c()

for(i in 1:length(d2)){
  up[i]<-quantile(d2[,i],0.975)
  lo[i]<-quantile(d2[,i],0.025)
  M[i]<-mean(d2[,i])
}

df<-data.frame(y=y$y, t=1:y$T, x=M, upper=up, lower=lo)
g<-ggplot(data = df, aes(x = t, y = y))
g<-g+geom_point()
g<-g+geom_errorbar(data = df, aes(ymin=lower, ymax=upper, x=t))
g<-g+geom_line(aes(x=t,y=x),data=df)
g

## トレンドモデル 二回差分のトレンド項を持つ、次のような状態空間モデルを考える。

tを時刻、\(x_t\)を状態、\(y_t\)を観測値として、 \[ y_t=x_t+\epsilon \] \[ x_t-x_{t-1}=x_{t-1}-x_{t-2}+\delta \] \[ \epsilon \sim N(0,\sigma_1) \] \[ \delta \sim N(0,\sigma_2) \]

推定するパラメータは状態\(x_t\)及び誤差の分散 \(\sigma_1,\sigma_2\)

data {
  int T; //サンプルサイズ
  real y[T]; //観測値
}

parameters {
  real x[T]; //状態
  real<lower=0> sigma_s; //状態のノイズの分散
  real<lower=0> sigma_o; //観測のノイズの分散
}

model {
  //観測方程式
  for (i in 1:T){
    y[i] ~ normal(x[i], sigma_o);
  }
  //状態方程式
  for (i in 3:T){
    x[i] ~ normal(2*x[i-1]-x[i-2], sigma_s);
  }
}

#ナイル川データ
Nile

## Time Series:
## Start = 1871 
## End = 1970 
## Frequency = 1 
##   [1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935 1110  994
##  [15] 1020  960 1180  799  958 1140 1100 1210 1150 1250 1260 1220 1030 1100
##  [29]  774  840  874  694  940  833  701  916  692 1020 1050  969  831  726
##  [43]  456  824  702 1120 1100  832  764  821  768  845  864  862  698  845
##  [57]  744  796 1040  759  781  865  845  944  984  897  822 1010  771  676
##  [71]  649  846  812  742  801 1040  860  874  848  890  744  749  838 1050
##  [85]  918  986  797  923  975  815 1020  906  901 1170  912  746  919  718
##  [99]  714  740

y<-list(y=as.numeric(Nile),T=length(Nile))
y

## $y
##   [1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935 1110  994
##  [15] 1020  960 1180  799  958 1140 1100 1210 1150 1250 1260 1220 1030 1100
##  [29]  774  840  874  694  940  833  701  916  692 1020 1050  969  831  726
##  [43]  456  824  702 1120 1100  832  764  821  768  845  864  862  698  845
##  [57]  744  796 1040  759  781  865  845  944  984  897  822 1010  771  676
##  [71]  649  846  812  742  801 1040  860  874  848  890  744  749  838 1050
##  [85]  918  986  797  923  975  815 1020  906  901 1170  912  746  919  718
##  [99]  714  740
## 
## $T
## [1] 100

library(rstan)
fit <- stan(file = 'trendmodel.stan', data = y,
                 iter = 4000, chains = 4)

## 
## SAMPLING FOR MODEL '8f4ae7eb030d6a862908f49755c46a85' NOW (CHAIN 1).
## 
## Gradient evaluation took 3.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.34 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 37.148 seconds (Warm-up)
##                36.1094 seconds (Sampling)
##                73.2574 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '8f4ae7eb030d6a862908f49755c46a85' NOW (CHAIN 2).
## 
## Gradient evaluation took 1.9e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 25.8401 seconds (Warm-up)
##                35.8858 seconds (Sampling)
##                61.726 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '8f4ae7eb030d6a862908f49755c46a85' NOW (CHAIN 3).
## 
## Gradient evaluation took 2.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 29.6418 seconds (Warm-up)
##                38.0186 seconds (Sampling)
##                67.6604 seconds (Total)
## 
## 
## SAMPLING FOR MODEL '8f4ae7eb030d6a862908f49755c46a85' NOW (CHAIN 4).
## 
## Gradient evaluation took 2.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 4000 [  0%]  (Warmup)
## Iteration:  400 / 4000 [ 10%]  (Warmup)
## Iteration:  800 / 4000 [ 20%]  (Warmup)
## Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Iteration: 4000 / 4000 [100%]  (Sampling)
## 
##  Elapsed Time: 27.8733 seconds (Warm-up)
##                28.7883 seconds (Sampling)
##                56.6616 seconds (Total)

## Warning: There were 72 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: There were 5996 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

#状態の事後分布
x<-rstan::extract(fit)$x
#状態の推定結果を可視化する
head(x)

##           
## iterations     [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
##       [1,] 1151.912 1138.675 1129.613 1116.637 1106.010 1095.848 1085.508
##       [2,] 1079.696 1096.497 1108.558 1129.019 1136.307 1138.898 1139.619
##       [3,] 1168.518 1161.621 1153.822 1147.106 1139.229 1130.365 1121.841
##       [4,] 1118.545 1120.045 1127.601 1131.744 1139.188 1150.168 1157.487
##       [5,] 1161.904 1163.828 1166.223 1164.673 1159.081 1149.857 1136.358
##       [6,] 1103.376 1093.929 1090.733 1086.864 1085.628 1085.923 1090.159
##           
## iterations     [,8]     [,9]    [,10]    [,11]    [,12]    [,13]    [,14]
##       [1,] 1079.646 1073.527 1068.425 1060.723 1053.614 1051.133 1048.388
##       [2,] 1139.012 1136.292 1127.919 1117.816 1106.167 1100.940 1089.598
##       [3,] 1114.284 1101.535 1088.570 1075.414 1065.065 1052.467 1041.608
##       [4,] 1164.360 1162.936 1164.005 1159.711 1155.027 1146.907 1137.121
##       [5,] 1125.806 1115.733 1104.280 1091.526 1084.913 1080.167 1073.720
##       [6,] 1093.446 1092.957 1090.809 1089.780 1087.912 1087.322 1080.002
##           
## iterations    [,15]    [,16]    [,17]    [,18]     [,19]     [,20]
##       [1,] 1052.193 1056.128 1057.198 1057.986 1061.8360 1063.5906
##       [2,] 1084.457 1079.219 1081.915 1093.302 1098.6079 1103.1756
##       [3,] 1034.306 1027.182 1017.666 1006.769  995.7331  985.7394
##       [4,] 1127.765 1115.848 1101.399 1084.494 1066.4675 1051.8163
##       [5,] 1066.279 1063.775 1059.457 1054.385 1046.2131 1037.6971
##       [6,] 1071.965 1066.021 1058.246 1050.538 1044.0669 1042.3748
##           
## iterations     [,21]     [,22]     [,23]     [,24]     [,25]     [,26]
##       [1,] 1060.5871 1057.2639 1053.2364 1042.6406 1032.6445 1017.5486
##       [2,] 1098.4349 1095.0868 1092.7199 1081.1253 1068.3897 1044.6050
##       [3,]  976.2464  966.0276  960.0635  954.3403  948.0269  939.6815
##       [4,] 1048.7982 1053.7165 1055.6097 1051.1702 1047.9714 1041.3988
##       [5,] 1033.7884 1029.9157 1027.3377 1023.2578 1013.0306 1001.6884
##       [6,] 1040.5707 1036.5167 1030.3623 1023.4962 1013.7124  998.0347
##           
## iterations     [,27]     [,28]    [,29]    [,30]    [,31]    [,32]
##       [1,]  998.5276  979.3911 963.8584 948.3173 936.7889 918.9724
##       [2,] 1024.4354  991.7257 967.7781 946.4270 917.5822 902.1012
##       [3,]  931.6988  926.7312 919.3501 911.0979 905.8389 897.6285
##       [4,] 1032.5858 1015.2703 999.2744 976.5877 952.8864 930.0904
##       [5,]  989.0876  975.5091 961.9451 947.5765 936.7332 927.8764
##       [6,]  985.3367  971.1221 958.6742 948.3490 937.9489 924.8826
##           
## iterations    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]    [,39]
##       [1,] 903.4111 883.5391 865.7180 848.0481 833.6703 824.1482 814.8831
##       [2,] 884.5276 864.9298 852.2904 848.2961 849.8624 849.0716 847.4615
##       [3,] 892.4018 890.0561 891.2037 890.6709 889.2986 888.8122 889.8567
##       [4,] 911.3095 899.1366 885.6676 878.7664 869.0291 859.0229 852.3670
##       [5,] 918.9502 909.8816 902.6129 895.0028 887.1163 878.5194 877.9684
##       [6,] 909.6029 891.7147 875.1743 859.5162 848.9953 838.9079 829.8826
##           
## iterations    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
##       [1,] 810.7338 807.2904 803.2117 799.5795 793.5651 791.9936 785.6289
##       [2,] 845.9335 845.5229 843.7040 835.2111 828.0422 826.7996 823.9406
##       [3,] 887.8791 885.0573 880.6889 878.0091 873.8922 872.0122 867.3729
##       [4,] 845.1886 839.9512 841.3411 837.0619 828.3861 817.1237 812.1087
##       [5,] 875.6876 872.8584 871.8226 867.6691 864.3318 859.0730 849.3860
##       [6,] 821.1311 816.4835 812.3594 806.7098 800.5131 797.0361 797.9078
##           
## iterations    [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]
##       [1,] 781.6003 778.9037 783.4383 787.7365 792.3042 799.4801 807.5065
##       [2,] 823.4446 822.9554 817.7881 814.8289 813.0558 810.1638 808.1586
##       [3,] 861.1998 857.0189 853.0120 851.1664 849.3152 846.6463 844.5387
##       [4,] 806.5194 794.0574 787.2644 779.4266 774.8014 765.4189 757.5524
##       [5,] 842.8422 838.8921 836.8531 831.8350 828.5366 828.1430 828.1316
##       [6,] 795.5870 794.1862 792.6495 793.6189 797.4821 802.3138 811.4986
##           
## iterations    [,54]    [,55]    [,56]    [,57]    [,58]    [,59]    [,60]
##       [1,] 815.8806 822.5602 829.5918 838.6114 842.0441 842.2426 838.0464
##       [2,] 812.2984 818.5740 824.3992 828.7304 836.7801 848.7685 862.7127
##       [3,] 843.9928 845.8369 848.2084 849.8476 851.3488 853.4175 851.4572
##       [4,] 751.5262 749.8805 751.2438 762.7057 762.4196 764.4859 775.9435
##       [5,] 830.3395 831.7906 834.4389 837.5315 839.8059 841.0317 843.0159
##       [6,] 818.7833 823.7254 829.3903 834.4325 840.3633 842.7943 846.3338
##           
## iterations    [,61]    [,62]    [,63]    [,64]    [,65]    [,66]    [,67]
##       [1,] 833.9893 826.8623 820.9975 814.0002 809.0675 806.4417 803.8717
##       [2,] 872.3049 876.7166 881.3386 881.4250 879.3655 885.6822 884.1773
##       [3,] 845.0239 836.7235 834.5200 832.6625 833.1105 832.8155 829.4165
##       [4,] 783.9937 792.3477 800.2023 813.1272 817.2055 809.8553 800.9833
##       [5,] 847.2243 852.4591 853.4193 856.2126 855.5755 853.6586 849.0259
##       [6,] 844.2030 841.5322 841.0181 841.3179 842.6814 848.1771 849.3083
##           
## iterations    [,68]    [,69]    [,70]    [,71]    [,72]    [,73]    [,74]
##       [1,] 796.9844 796.7716 796.6779 798.9893 804.0298 812.9963 830.7010
##       [2,] 878.0300 868.8934 852.0444 833.0466 823.1163 818.7826 812.3427
##       [3,] 826.3337 828.5917 833.5832 836.2110 838.9027 842.7931 845.3681
##       [4,] 795.4023 794.2715 796.8597 803.4575 807.2521 812.4845 822.9450
##       [5,] 842.0500 837.4591 835.6377 833.5402 834.4540 832.4625 832.7407
##       [6,] 850.7037 853.3505 853.7928 853.1767 856.6335 857.1905 854.7183
##           
## iterations    [,75]    [,76]    [,77]    [,78]    [,79]    [,80]    [,81]
##       [1,] 845.0724 855.6097 866.2605 881.0214 892.9593 902.3857 909.8784
##       [2,] 816.7183 812.6026 815.8092 814.0597 806.6927 805.3715 808.0404
##       [3,] 848.7692 853.2094 855.6922 856.0358 857.2948 860.1728 864.6937
##       [4,] 821.1742 818.5064 835.2823 844.8418 850.6160 854.6553 857.5068
##       [5,] 835.1573 839.3570 840.2635 842.1349 844.2863 846.9384 850.2459
##       [6,] 850.3042 845.6502 845.6504 845.5643 849.4115 849.0322 848.9174
##           
## iterations    [,82]    [,83]    [,84]    [,85]    [,86]    [,87]    [,88]
##       [1,] 923.0105 933.0436 944.4725 953.8134 958.7521 960.6812 959.9637
##       [2,] 814.2269 823.5158 828.9985 829.2456 834.4347 846.6707 860.6411
##       [3,] 871.0863 875.4631 878.5503 881.6220 883.9119 886.2900 889.1937
##       [4,] 863.9558 877.4412 890.5246 896.5218 900.9259 895.9126 891.9759
##       [5,] 855.6405 859.5346 865.8680 873.7774 884.0166 890.6127 893.4777
##       [6,] 848.5081 847.8315 849.9815 848.6301 847.0390 841.8645 838.1803
##           
## iterations    [,89]    [,90]    [,91]    [,92]    [,93]    [,94]    [,95]
##       [1,] 962.6932 961.5718 956.2538 947.1501 930.8453 908.4422 878.9366
##       [2,] 869.0474 889.4999 897.1590 891.8792 881.0941 863.4229 845.3141
##       [3,] 890.2752 894.4677 897.3549 899.0677 902.6687 904.1823 906.4412
##       [4,] 897.4825 883.5493 878.1161 873.4262 857.3368 838.5133 820.4832
##       [5,] 894.0014 893.7118 890.4629 883.0152 873.1217 858.9095 845.9962
##       [6,] 840.6796 839.7442 839.5869 840.8020 839.4240 841.9605 845.4881
##           
## iterations    [,96]    [,97]    [,98]    [,99]   [,100]
##       [1,] 844.5893 815.0946 786.7177 753.4059 722.6501
##       [2,] 825.4614 798.5900 768.5616 739.0253 717.9721
##       [3,] 907.8436 906.1885 902.3339 897.1697 892.2848
##       [4,] 799.9522 787.7580 772.9346 757.0092 733.1642
##       [5,] 833.3056 820.0104 806.1767 795.6975 787.1935
##       [6,] 845.4052 843.4260 841.0926 829.8172 817.6404

d2<-data.frame(x)
up<-c()
lo<-c()
M<-c()

for(i in 1:length(d2)){
  up[i]<-quantile(d2[,i],0.975)
  lo[i]<-quantile(d2[,i],0.025)
  M[i]<-mean(d2[,i])
}

df<-data.frame(y=y$y, t=1:y$T, x=M, upper=up, lower=lo)
g<-ggplot(data = df, aes(x = t, y = y))
g<-g+geom_point()
g<-g+geom_errorbar(data = df, aes(ymin=lower, ymax=upper, x=t))
g<-g+geom_line(aes(x=t,y=x),data=df)
g